{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Important details of the first order formulation\n",
    "This notebook here is to test the first order ODE solver\n",
    "\n",
    "The main thing to worry about is the fact that we need to do some mode sorting in order to properly evaluate the eigenmodes for the $E$ and $H$ fields (so it matches with the 2nd order formulation, which automatically does 'mode sorting').\n",
    "\n",
    "Beyond that, we should be able to extract both the $W$ and $V$ eigenmodes for the H and E field eigenmodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\zhaon\\Anaconda3\\lib\\site-packages\\numpy\\core\\numeric.py:501: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## same as the analytic case but with the fft\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.linalg import cond\n",
    "import cmath;\n",
    "from scipy.fftpack import fft, fftfreq, fftshift, rfft\n",
    "from scipy.fftpack import dst, idst\n",
    "from scipy.linalg import expm\n",
    "from scipy import linalg as LA\n",
    "import random\n",
    "\n",
    "# Moharam et. al Formulation for stable and efficient implementation for RCWA\n",
    "plt.close(\"all\")\n",
    "'''\n",
    "1D TE implementation of PLANAR DIFFRACTiON...the easy case\n",
    "only: sign convention is exp(-ikr) (is the positive propagating wave), so loss is +  not - \n",
    "source for fourier decomps is from the paper: Formulation for stable and efficient implementation of\n",
    "the rigorous coupled-wave analysis of binary gratings by Moharam et. al\n",
    "'''\n",
    "np.set_printoptions(precision = 4)\n",
    "\n",
    "def grating_fourier_harmonics(order, fill_factor, n_ridge, n_groove):\n",
    "    \"\"\" function comes from analytic solution of a step function in a finite unit cell\"\"\"\n",
    "    #n_ridge = index of refraction of ridge (should be dielectric)\n",
    "    #n_ridge = index of refraction of groove (air)\n",
    "    #n_ridge has fill_factor\n",
    "    #n_groove has (1-fill_factor)\n",
    "    # there is no lattice constant here, so it implicitly assumes that the lattice constant is 1...which is not good\n",
    "\n",
    "    if(order == 0):\n",
    "        return n_ridge**2*fill_factor + n_groove**2*(1-fill_factor);\n",
    "    else:\n",
    "        #should it be 1-fill_factor or fill_factor?, should be fill_factor\n",
    "        return(n_ridge**2 - n_groove**2)*np.sin(np.pi*order*(fill_factor))/(np.pi*order);\n",
    "\n",
    "def grating_fourier_array(num_ord, fill_factor, n_ridge, n_groove):\n",
    "    \"\"\" what is a convolution in 1D \"\"\"\n",
    "    fourier_comps = list();\n",
    "    for i in range(-num_ord, num_ord+1):\n",
    "        fourier_comps.append(grating_fourier_harmonics(i, fill_factor, n_ridge, n_groove));\n",
    "    return fourier_comps;\n",
    "\n",
    "def fourier_reconstruction(x, period, num_ord, n_ridge, n_groove, fill_factor = 0.5):\n",
    "    index = np.arange(-num_ord, num_ord+1);\n",
    "    f = 0;\n",
    "    for n in index:\n",
    "        coef = grating_fourier_harmonics(n, fill_factor, n_ridge, n_groove);\n",
    "        f+= coef*np.exp(cmath.sqrt(-1)*np.pi*n*x/period);\n",
    "        #f+=coef*np.cos(np.pi*n*x/period)\n",
    "    return f;\n",
    "\n",
    "def fourier_reconstruction_general(x, period, num_ord, coefs):\n",
    "    '''\n",
    "    overloading odesn't work in python...fun fact, since it is dynamically typed (vs statically typed)\n",
    "    :param x:\n",
    "    :param period:\n",
    "    :param num_ord:\n",
    "    :param coefs:\n",
    "    :return:\n",
    "    '''\n",
    "    index = np.arange(-num_ord, num_ord+1);\n",
    "    f = 0; center = int(len(coefs)/2); #no offset\n",
    "    for n in index:\n",
    "        coef = coefs[center+n];\n",
    "        f+= coef*np.exp(cmath.sqrt(-1)*2*np.pi*n*x/period);\n",
    "    return f;\n",
    "\n",
    "def grating_fft(eps_r):\n",
    "    assert len(eps_r.shape) == 2\n",
    "    assert eps_r.shape[1] == 1;\n",
    "    #eps_r: discrete 1D grid of the epsilon profile of the structure\n",
    "    fourier_comp = np.fft.fftshift(np.fft.fft(eps_r, axis = 0)/eps_r.shape[0]);\n",
    "    #ortho norm in fft will do a 1/sqrt(n) scaling\n",
    "    return np.squeeze(fourier_comp);\n",
    "\n",
    "# plt.plot(x, np.real(fourier_reconstruction(x, period, 1000, 1,np.sqrt(12), fill_factor = 0.1)));\n",
    "# plt.title('check that the analytic fourier series works')\n",
    "# #'note that the lattice constant tells you the length of the ridge'\n",
    "# plt.show()\n",
    "\n",
    "L0 = 1e-6;\n",
    "e0 = 8.854e-12;\n",
    "mu0 = 4*np.pi*1e-8;\n",
    "fill_factor = 0.3; # 50% of the unit cell is the ridge material\n",
    "\n",
    "\n",
    "num_ord = 10; #INCREASING NUMBER OF ORDERS SEEMS TO CAUSE THIS THING TO FAIL, to many orders induce evanescence...particularly\n",
    "               # when there is a small fill factor\n",
    "PQ = 2*num_ord+1;\n",
    "indices = np.arange(-num_ord, num_ord+1)\n",
    "\n",
    "n_ridge = 3.48; #1-3.48j; #3.48;              # ridge\n",
    "n_groove =  1;                # groove (unit-less)\n",
    "lattice_constant = 0.7;  # SI units\n",
    "# we need to be careful about what lattice constant means\n",
    "# in the gaylord paper, lattice constant exactly means (0, L) is one unit cell\n",
    "\n",
    "\n",
    "d = 0.46;               # thickness, SI units\n",
    "Nx = 2*256;\n",
    "eps_r = n_groove**2*np.ones((2*Nx, 1)); #put in a lot of points in eps_r\n",
    "border = int(2*Nx*fill_factor);\n",
    "eps_r[0:border] = n_ridge**2;\n",
    "fft_fourier_array = grating_fft(eps_r);\n",
    "x = np.linspace(-lattice_constant,lattice_constant,1000);\n",
    "period = lattice_constant;\n",
    "fft_reconstruct = fourier_reconstruction_general(x, period, num_ord, fft_fourier_array);\n",
    "\n",
    "fourier_array_analytic = grating_fourier_array(Nx, fill_factor, n_ridge, n_groove);\n",
    "analytic_reconstruct = fourier_reconstruction(x, period, num_ord, n_ridge, n_groove, fill_factor)\n",
    "\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(np.real(fft_fourier_array[Nx-20:Nx+20]), linewidth=2)\n",
    "plt.plot(np.real(fourier_array_analytic[Nx-20:Nx+20]));\n",
    "plt.legend(('fft', 'analytic'))\n",
    "plt.show()\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(x,fft_reconstruct)\n",
    "plt.plot(x,analytic_reconstruct);\n",
    "plt.legend(['fft', 'analytic'])\n",
    "plt.show()\n",
    "\n",
    "## simulation parameters\n",
    "theta = (0)*np.pi/180;\n",
    "\n",
    "\n",
    "## construct permittivity harmonic components E\n",
    "#fill factor = 0 is complete dielectric, 1 is air\n",
    "\n",
    "##construct convolution matrix\n",
    "E = np.zeros((2 * num_ord + 1, 2 * num_ord + 1)); E = E.astype('complex')\n",
    "p0 = Nx; #int(Nx/2);\n",
    "p_index = np.arange(-num_ord, num_ord + 1);\n",
    "q_index = np.arange(-num_ord, num_ord + 1);\n",
    "fourier_array = fft_fourier_array;#fourier_array_analytic;\n",
    "detected_pffts = np.zeros_like(E);\n",
    "for prow in range(2 * num_ord + 1):\n",
    "    # first term locates z plane, 2nd locates y coumn, prow locates x\n",
    "    row_index = p_index[prow];\n",
    "    for pcol in range(2 * num_ord + 1):\n",
    "        pfft = p_index[prow] - p_index[pcol];\n",
    "        detected_pffts[prow, pcol] = pfft;\n",
    "        E[prow, pcol] = fourier_array[p0 + pfft];  # fill conv matrix from top left to top right\n",
    "\n",
    "## IMPORTANT TO NOTE: the indices for everything beyond this points are indexed from -num_ord to num_ord+1\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Calculating the Poynting Vector of a Sum of Plane Waves\n",
    "Each plane wave has a contributing amplitude. Do we sum all the plane waves then calculate the Poynting vector or can we sum the Poynting vectors of the individual waves? The issue is that a product of a sum is not the same as the sum of products. Presently, I've tried both. It only appears to work in the simple case of 0 orders (or just plane waves).\n",
    "\n",
    "Technically, each plane has has a different propagation in $k_z$ indexed by their Fourier index, which means we cna't deal with the components by themselves, we'd have to add in their phase. We could also just set $z=0$ to avoid this nasty problem.\n",
    "\n",
    "Secret Sauce: Obviously with the Fourier order decomposition in $x$, some of the fourier $k_{zi}$ components may be evanescent. These components should have no contribution to the Poynting vector (especially in the far-field).\n",
    "\n",
    "### some debugging notes\n",
    "the long wavelength case should be error-free because there is little scattering into higher diffraction orders... but it appears to be the opposite in our test case below. Mathematically, however, we can see it in the fact that we have $m_i \\lambda$ when we determine our kx_array, which means larger wavelengths have larger contributions at higher orders, which is a weird sounding statement.\n",
    "\n",
    "In the 2nd order, stable formulation, we simply have REAL eigenvalues (if the system is Hermitian). When we square root this, if we have a purely negative number, then doing the square root should still only yield a single sign in the imaginary part as well.\n",
    "\n",
    "It is clear that solving the first order and second order eigenvalue problems lead to the same set of eigenvalues, but the choice selected by the second order problem seems weird. Specifically, while the real part is always consistently one sign (which means sorting would work), the imaginary part selection is not all the same sign, particularly the eigenvalues with effectively zero real $k$. The problem with this has been determined to be a numerical artifact (see below). One proposed solution is to zero out the real part if it's below numerical precision (but this feels very unsatisfactory)\n",
    "\n",
    "Empirically, it seems that sorting the larger block eigenvalues FAILS period.\n",
    "\n",
    "In our current code implementation, we have the X matrix as np.exp(-k_0 *Q *d). So we neg out everything. However, the fact that some of the imaginary parts flip sign doesn't mean this part is wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## plot snippets\n",
    "## plot verifies that all eigenvalues extracted by beigenvals matches with one in Q\n",
    "#     plt.subplot(121)\n",
    "#     for b in beigenvals:\n",
    "#         plt.axhline(np.imag(b))\n",
    "#     plt.plot(np.imag(np.diag(Q)),'.g-', markersize = 10)\n",
    "#     plt.title(str((len(beigenvals), len(np.diag(Q)))))\n",
    "#     plt.subplot(122)\n",
    "#     for b in beigenvals:\n",
    "#         plt.axhline(np.real(b))\n",
    "#     plt.plot(np.real(np.diag(Q)),'.g-', markersize = 10)\n",
    "#     plt.title('real')\n",
    "#     plt.show()\n",
    "\n",
    "## snippet to check maching square of eigenvals\n",
    "#     for b in np.abs(np.square(beigenvals)):\n",
    "# #         plt.axhline(b);\n",
    "\n",
    "# #     plt.plot(sorted(np.abs(eigenvals)));\n",
    "# #     plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-3-aadf50ac8741>, line 130)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-3-aadf50ac8741>\"\u001b[1;36m, line \u001b[1;32m130\u001b[0m\n\u001b[1;33m    V = W@Q; #H modes\u001b[0m\n\u001b[1;37m         ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "## alternate construction of 1D convolution matrix\n",
    "spectra = list();\n",
    "spectra_T = list();\n",
    "I = np.identity(2 * num_ord + 1)\n",
    "PQ = 2*num_ord+1;\n",
    "zeros = np.zeros((PQ, PQ))\n",
    "# E is now the convolution of fourier amplitudes\n",
    "wavelength_scan = np.linspace(0.5,2.3,300)\n",
    "\n",
    "for wvlen in wavelength_scan:\n",
    "    j = cmath.sqrt(-1);\n",
    "    lam0 = wvlen;     k0 = 2 * np.pi / lam0; #free space wavelength in SI units\n",
    "    print('wavelength: ' + str(wvlen));\n",
    "    ## =====================STRUCTURE======================##\n",
    "\n",
    "    ## Region I: reflected region (half space)\n",
    "    n1 = 1;#cmath.sqrt(-1)*1e-12; #apparently small complex perturbations are bad in Region 1, these shouldn't be necessary\n",
    "\n",
    "    ## Region 2; transmitted region\n",
    "    n2 = 1;\n",
    "\n",
    "    #from the kx_components given the indices and wvln\n",
    "    kx_array = k0*(n1*np.sin(theta) + indices*(lam0 / lattice_constant)); #0 is one of them, k0*lam0 = 2*pi\n",
    "    k_xi = kx_array;\n",
    "    ## IMPLEMENT SCALING: these are the fourier orders of the x-direction decomposition.\n",
    "    KX = np.diag(kx_array/k0);\n",
    "    KX2 = np.diag(np.power((k_xi/k0),2)); #singular since we have a n=0, m= 0 order and incidence is normal\n",
    "    #print(KX2)\n",
    "    KZ2 = ((n1)**2 - np.diag(KX2)).astype('complex');\n",
    "    KZ = np.sqrt(KZ2)\n",
    "#     print('KZ2: '+str(KZ2)+': imag Kz: '+str(KZ))\n",
    "#     KZ_mask = np.imag(KZ)<1e-2;\n",
    "#     KZ_mask = np.expand_dims(KZ_mask, axis = 1)\n",
    "#     print(KZ_mask.shape)\n",
    "    \n",
    "    ## construct matrix of Gamma^2 ('constant' term in ODE):\n",
    "    A = KX2 - E; #conditioning of this matrix is not bad, A SHOULD BE SYMMETRIC\n",
    "    #sum of a symmetric matrix and a diagonal matrix should be symmetric;\n",
    "    AO = np.block([[zeros, I],[A, zeros]])\n",
    "    beigenvals, bigW = LA.eig(AO);\n",
    "    #print('conditioning of big block: '+str(np.linalg.cond(AO)))\n",
    "    ## SORTING IS REQUIRED\n",
    "    #sorting procedure is to order by smallest imaginary part to largest imaginary part\n",
    "    # We actually have to calculate poynting vector (the modes are Ey, Hx) and Sz = -Ey*Hx\n",
    "    # since we are dealing with a normalized H, we need to divide by 1/j;\n",
    "    Ey_modes = bigW[0:PQ, :];\n",
    "    Hx_modes = bigW[PQ:, :]\n",
    "    \n",
    "\n",
    "    EyAmp = np.sum(Ey_modes, axis = 0); #amplitude\n",
    "    HxAmp = np.sum(Hx_modes, axis = 0);\n",
    "    #EyAmp = np.sum(Ey_modes*KZ_mask, axis = 0); #amplitude\n",
    "    #HxAmp = np.sum(Hx_modes*KZ_mask, axis = 0);\n",
    "    \n",
    "    ## right now it works for a uniform slab, any orders...which means it \n",
    "    Sz = (1/j)*EyAmp*(HxAmp);\n",
    "#     print(Sz.shape)\n",
    "    sorted_indices_poynting = np.argsort(Sz)\n",
    "    sorted_indices = np.argsort(np.real(beigenvals))\n",
    "    #print(Sz[sorted_indices])\n",
    "    \n",
    "    sorted_eigenmodes = bigW[:,sorted_indices ];\n",
    "    #print(sorted_eigenmodes)\n",
    "    #adding real and imaginary parts seems to work...\n",
    "    sorted_eigenvals = beigenvals[sorted_indices]\n",
    "    Wp = sorted_eigenmodes[0:PQ:,0:PQ]\n",
    "    eigenvals_wp = (sorted_eigenvals[0:PQ]);\n",
    "    #print('sorted beigenvals: '+str(sorted_eigenvals))\n",
    "    #print('extracted eigenvals: '+str(eigenvals_wp))\n",
    "\n",
    "    #print(eigenvals_wp.shape)\n",
    "    #plt.plot(np.imag(sorted_eigenvals));\n",
    "    #plt.show()\n",
    "    ##\n",
    "    # when we calculate eigenvals, how do we know the eigenvals correspond to each particular fourier order?\n",
    "    eigenvals, W = LA.eig(A); #A should be symmetric or hermitian\n",
    "    #print('sorted eigenvals: '+str(sorted(np.sqrt(eigenvals))))\n",
    "\n",
    "#     plt.subplot(121);\n",
    "#     plt.imshow(np.abs(W));\n",
    "#     plt.subplot(122);\n",
    "#     plt.imshow(np.abs(bigW));\n",
    "#     plt.show()\n",
    "    \n",
    "    #W = Wp; \n",
    "    #V = sorted_eigenmodes[0:PQ, :]\n",
    "    #we should be gauranteed that all eigenvals are REAL\n",
    "    eigenvals = eigenvals.astype('complex');\n",
    "    \n",
    "    Q = np.diag(np.sqrt(eigenvals)); #Q should only be positive square root of eigenvals\n",
    "    #real parts match, but the imaginaries don't\n",
    "    #Q = np.diag(eigenvals_wp);\n",
    "    #plt.plot(sorted(np.abs(beigenvals)),'.-');\n",
    "#     colors = [(random.random(),random.random(),random.random()) for i in range(len(beigenvals))]\n",
    "#     plt.subplot(131)\n",
    "#     plt.plot(np.real(beigenvals), np.imag(beigenvals), '.'); plt.title('1st')\n",
    "#     plt.ylim([-4,4])\n",
    "#     plt.subplot(132)\n",
    "#     plt.plot(np.real(np.diag(Qn)), np.imag(np.diag(Qn)), '.'); plt.title('2nd')\n",
    "#     plt.ylim([-4,4])\n",
    "#     plt.subplot(133);\n",
    "#     plt.plot(np.real(np.sqrt(eigenvals)), np.imag(np.sqrt(eigenvals)), '.g'); plt.title('squared 2nd')\n",
    "#     plt.ylim([-4,4])\n",
    "\n",
    "    plt.show();\n",
    "#     plt.subplot(121)\n",
    "#     for b in beigenvals:\n",
    "#         plt.axhline(np.imag(b), color='r')\n",
    "#     for b in np.sqrt(eigenvals):\n",
    "#         plt.axhline(np.imag(b))\n",
    "\n",
    "#     plt.plot((np.imag(np.diag(Q))),'.g-', markersize = 10)\n",
    "#     plt.plot((np.imag(np.diag(Qn))), '.c-', markersize = 10)\n",
    "#     plt.title(str((len(beigenvals), len(np.diag(Q)))))\n",
    "    \n",
    "#     plt.subplot(122)\n",
    "    \n",
    "#     for b in beigenvals:\n",
    "#         plt.axhline(np.real(b), color = 'r')\n",
    "#     for b in np.sqrt(eigenvals):\n",
    "#         plt.axhline(np.real(b))\n",
    "        \n",
    "#     plt.plot((np.real(np.diag(Q))),'.g-', markersize = 10)\n",
    "#     plt.title('real')\n",
    "\n",
    "#     plt.plot((np.real(np.diag(Qn))), '.c-')\n",
    "#     plt.show()\n",
    "\n",
    "\n",
    "    V = W@Q; #H modes\n",
    "\n",
    "    ## THIS ENFORCES a PARTICULAR SIGN...\n",
    "    X = np.diag(np.exp(-k0*np.diag(Q)*d)); #this is poorly conditioned because exponentiation\n",
    "    ## pointwise exponentiation vs exponentiating a matrix\n",
    "\n",
    "    ## observation: almost everything beyond this point is worse conditioned\n",
    "    k_I = k0**2*(n1**2 - (k_xi/k0)**2);                 #k_z in reflected region k_I,zi\n",
    "    k_II = k0**2*(n2**2 - (k_xi/k0)**2);   #k_z in transmitted region\n",
    "    k_I = k_I.astype('complex'); k_I = np.sqrt(k_I);\n",
    "    k_II = k_II.astype('complex'); k_II = np.sqrt(k_II);\n",
    "    Y_I = np.diag(k_I/k0);\n",
    "    Y_II = np.diag(k_II/k0);\n",
    "\n",
    "\n",
    "    delta_i0 = np.zeros((len(kx_array),1));\n",
    "    delta_i0[num_ord] = 1;\n",
    "    n_delta_i0 = delta_i0*j*n1*np.cos(theta); #this is a VECTOR\n",
    "\n",
    "    ## design auxiliary variables: SEE derivation in notebooks: RCWA_note.ipynb\n",
    "    # we want to design the computation to avoid operating with X, particularly with inverses\n",
    "    # since X is the worst conditioned thing\n",
    "\n",
    "    Wi = np.linalg.inv(W);\n",
    "    Vi = np.linalg.inv(V);\n",
    "    Oi = 0.5*np.block([[Wi, Vi],[Wi, -Vi]])\n",
    "    f = I;\n",
    "    g = j*Y_II; #all matrices\n",
    "    fg = np.concatenate((f,g),axis = 0)\n",
    "    #ab = np.matmul(np.linalg.inv(O),fg);\n",
    "    # ab = np.matmul(Oi, fg);\n",
    "    # a = ab[0:PQ,:];\n",
    "    # b = ab[PQ:,:];\n",
    "\n",
    "    a = 0.5*(Wi+j*Vi@Y_II);\n",
    "    b = 0.5*(Wi-j*Vi@Y_II);\n",
    "    fbiX = np.matmul(np.linalg.inv(b),X)\n",
    "\n",
    "    #altTerm = (a@X@X@b); #not well conditioned and I-altTermis is also poorly conditioned.\n",
    "    #print(np.linalg.cond(I-np.linalg.inv(altTerm)))\n",
    "    #print(np.linalg.cond(X@b)); #not well conditioned.\n",
    "\n",
    "    term = X@a@fbiX; # THIS IS SHITTILY CONDITIONED\n",
    "    # print((np.linalg.cond(X), np.linalg.cond(term)))\n",
    "    # print(np.linalg.cond(I+term)); #but this is EXTREMELY WELL CONDITIONED.\n",
    "    f = np.matmul(W, I+term);\n",
    "    g = np.matmul(V,-I+term);\n",
    "    T = np.linalg.inv(j*np.matmul(Y_I,f)+g);\n",
    "    T = np.matmul(T,(np.matmul(j*Y_I,delta_i0)+n_delta_i0));\n",
    "    R = np.matmul(f,T)-delta_i0;\n",
    "    T = np.matmul(fbiX, T)\n",
    "\n",
    "    ## calculate diffraction efficiencies\n",
    "    #I would expect this number to be real...\n",
    "    DE_ri = R*np.conj(R)*np.real(np.expand_dims(k_I,1))/(k0*n1*np.cos(theta));\n",
    "    DE_ti = T*np.conj(T)*np.real(np.expand_dims(k_II,1))/(k0*n1*np.cos(theta));\n",
    "    #print(np.sum(DE_ri))\n",
    "\n",
    "    #print(np.sum(DE_ri))\n",
    "    spectra.append(np.sum(DE_ri)); #spectra_T.append(T);\n",
    "    spectra_T.append(np.sum(DE_ti))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'wavelength_scan' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-4-05ca4470c35b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwavelength_scan\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mspectra\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwavelength_scan\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mspectra_T\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'wavelength (microns)'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'wavelength_scan' is not defined"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.figure();\n",
    "plt.plot(wavelength_scan, np.abs(spectra));\n",
    "plt.plot(wavelength_scan, np.abs(spectra_T))\n",
    "plt.xlabel('wavelength (microns)')\n",
    "plt.ylabel('R/T')\n",
    "plt.title('sample RCWA spectra for a 1D grating')\n",
    "plt.legend(['reflection', 'transmission'])\n",
    "# plt.axhline(((3.48-1)/(3.48+1))**2,xmin=0, xmax = max(wavelength_scan))\n",
    "# plt.axhline(((3.48-1)/(3.48+1)),xmin=0, xmax = max(wavelength_scan), color='r')\n",
    "plt.ylim([0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## comparison to the 2nd order, stable formulation\n",
    "# also the first real test of the function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Curious Numerical Observation\n",
    "Something which plagues the first order formulation is that if we take a number,say -5. Taking the square root would give $-\\sqrt{5}$. however, if we add in a tiny COMPLEX perturbation of negative sign, we can completely flip the sign. This must be a numerical artifact of numerical imprecision that happens in the stable formulation when we take the square root, but WHY is the 2nd order formulation still stable?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2.2360679774997896e-17-2.23606797749979j)\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "a = -5 - j*1e-16;\n",
    "print(cmath.sqrt(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ONLY FIRST ORDER RAW CODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wavelength: 0.5\n",
      "conditioning of big block: 113.0809969549768\n",
      "[0, 1, 4, 5, 8, 9, 12, 13, 14, 18, 22, 23, 26, 30, 31, 32, 33, 36, 37, 40, 41]\n",
      "[ 2  3  6  7 10 11 15 16 17 19 20 21 24 25 27 28 29 34 35 38 39]\n",
      "(21, 21)\n",
      "wavelength: 0.5947368421052631\n",
      "conditioning of big block: 314.1987638798421\n",
      "[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 18, 21, 23, 24, 28, 29, 32, 33, 34, 38, 39]\n",
      "[ 2  3  6  7 10 11 14 15 19 20 22 25 26 27 30 31 35 36 37 40 41]\n",
      "(21, 21)\n",
      "wavelength: 0.6894736842105262\n",
      "conditioning of big block: 326.64766202236063\n",
      "[2, 3, 6, 7, 10, 11, 14, 15, 21, 22, 23, 24, 25, 26, 30, 31, 34, 35, 36, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 16 17 18 19 20 27 28 29 32 33 37 38 39]\n",
      "(21, 21)\n",
      "wavelength: 0.7842105263157895\n",
      "conditioning of big block: 175.54605652984324\n",
      "[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 18, 19, 25, 29, 32, 33, 34, 35, 37, 38, 39]\n",
      "[ 2  3  6  7 10 11 14 15 20 21 22 23 24 26 27 28 30 31 36 40 41]\n",
      "(21, 21)\n",
      "wavelength: 0.8789473684210526\n",
      "conditioning of big block: 319.51992816229966\n",
      "[2, 3, 6, 7, 10, 11, 14, 15, 18, 19, 22, 23, 24, 28, 29, 30, 34, 35, 38, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 16 17 20 21 25 26 27 31 32 33 36 37 39]\n",
      "(21, 21)\n",
      "wavelength: 0.9736842105263157\n",
      "conditioning of big block: 676.0122836735518\n",
      "[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 30, 31, 34, 38, 39, 40]\n",
      "[ 2  3  6  7 10 11 14 15 18 19 22 23 26 27 29 32 33 35 36 37 41]\n",
      "(21, 21)\n",
      "wavelength: 1.068421052631579\n",
      "conditioning of big block: 3540.421923989473\n",
      "[0, 1, 2, 3, 6, 7, 10, 11, 14, 15, 18, 19, 24, 25, 28, 30, 31, 32, 33, 34, 40]\n",
      "[ 4  5  8  9 12 13 16 17 20 21 22 23 26 27 29 35 36 37 38 39 41]\n",
      "(21, 21)\n",
      "wavelength: 1.163157894736842\n",
      "conditioning of big block: 1552.8877564970778\n",
      "[2, 3, 6, 7, 10, 11, 14, 15, 18, 19, 22, 23, 26, 27, 29, 30, 31, 32, 33, 38, 40]\n",
      "[ 0  1  4  5  8  9 12 13 16 17 20 21 24 25 28 34 35 36 37 39 41]\n",
      "(21, 21)\n",
      "wavelength: 1.2578947368421052\n",
      "conditioning of big block: 717.0211857895337\n",
      "[0, 1, 4, 5, 8, 9, 12, 13, 14, 15, 18, 19, 24, 25, 28, 30, 31, 32, 33, 34, 36]\n",
      "[ 2  3  6  7 10 11 16 17 20 21 22 23 26 27 29 35 37 38 39 40 41]\n",
      "(21, 21)\n",
      "wavelength: 1.3526315789473684\n",
      "conditioning of big block: 11125.338847198258\n",
      "[2, 3, 6, 7, 10, 11, 14, 15, 18, 19, 22, 23, 26, 27, 28, 32, 33, 36, 37, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 16 17 20 21 24 25 29 30 31 34 35 38 39]\n",
      "(21, 21)\n",
      "wavelength: 1.4473684210526314\n",
      "conditioning of big block: 570.8495898819775\n",
      "[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 29, 30, 31, 32, 33, 36, 37]\n",
      "[ 2  3  6  7 10 11 14 15 18 19 22 23 26 27 28 34 35 38 39 40 41]\n",
      "(21, 21)\n",
      "wavelength: 1.5421052631578946\n",
      "conditioning of big block: 481.1138035671346\n",
      "[2, 3, 6, 7, 10, 11, 14, 15, 18, 19, 24, 25, 26, 27, 32, 34, 35, 36, 37, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 16 17 20 21 22 23 28 29 30 31 33 38 39]\n",
      "(21, 21)\n",
      "wavelength: 1.6368421052631579\n",
      "conditioning of big block: 542.5622063132746\n",
      "[2, 3, 6, 7, 10, 11, 16, 17, 20, 21, 24, 25, 26, 27, 33, 34, 35, 36, 37, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 14 15 18 19 22 23 28 29 30 31 32 38 39]\n",
      "(21, 21)\n",
      "wavelength: 1.7315789473684209\n",
      "conditioning of big block: 607.6760071265111\n",
      "[0, 1, 2, 3, 6, 7, 12, 13, 14, 15, 18, 19, 22, 23, 28, 29, 30, 31, 33, 34, 35]\n",
      "[ 4  5  8  9 10 11 16 17 20 21 24 25 26 27 32 36 37 38 39 40 41]\n",
      "(21, 21)\n",
      "wavelength: 1.8263157894736841\n",
      "conditioning of big block: 676.4547935816985\n",
      "[2, 3, 6, 7, 10, 11, 16, 17, 20, 21, 24, 25, 26, 27, 32, 36, 37, 38, 39, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 14 15 18 19 22 23 28 29 30 31 33 34 35]\n",
      "(21, 21)\n",
      "wavelength: 1.9210526315789471\n",
      "conditioning of big block: 748.8982489308878\n",
      "[2, 3, 6, 7, 10, 11, 16, 17, 20, 21, 24, 25, 26, 27, 32, 36, 37, 38, 39, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 14 15 18 19 22 23 28 29 30 31 33 34 35]\n",
      "(21, 21)\n",
      "wavelength: 2.0157894736842104\n",
      "conditioning of big block: 825.0061268198372\n",
      "[0, 1, 4, 5, 8, 9, 10, 11, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31, 32, 34, 35]\n",
      "[ 2  3  6  7 12 13 16 17 18 19 24 25 26 27 33 36 37 38 39 40 41]\n",
      "(21, 21)\n",
      "wavelength: 2.110526315789474\n",
      "conditioning of big block: 904.7782334191904\n",
      "[2, 3, 6, 7, 10, 11, 16, 17, 18, 19, 24, 25, 26, 27, 32, 36, 37, 38, 39, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 14 15 20 21 22 23 28 29 30 31 33 34 35]\n",
      "(21, 21)\n",
      "wavelength: 2.205263157894737\n",
      "conditioning of big block: 988.2144146022588\n",
      "[2, 3, 6, 7, 10, 11, 16, 17, 18, 19, 24, 25, 26, 27, 33, 36, 37, 38, 39, 40, 41]\n",
      "[ 0  1  4  5  8  9 12 13 14 15 20 21 22 23 28 29 30 31 32 34 35]\n",
      "(21, 21)\n",
      "wavelength: 2.3\n",
      "conditioning of big block: 1075.3145466133653\n",
      "[0, 1, 4, 5, 8, 9, 10, 11, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31, 32, 33, 34]\n",
      "[ 2  3  6  7 12 13 16 17 18 19 24 25 26 27 35 36 37 38 39 40 41]\n",
      "(21, 21)\n"
     ]
    }
   ],
   "source": [
    "## alternate construction of 1D convolution matrix\n",
    "spectra = list();\n",
    "spectra_T = list();\n",
    "I = np.identity(2 * num_ord + 1)\n",
    "PQ = 2*num_ord+1;\n",
    "zeros = np.zeros((PQ, PQ))\n",
    "# E is now the convolution of fourier amplitudes\n",
    "wavelength_scan = np.linspace(0.5,2.3,20)\n",
    "\n",
    "for wvlen in wavelength_scan:\n",
    "    j = cmath.sqrt(-1);\n",
    "    lam0 = wvlen;     k0 = 2 * np.pi / lam0; #free space wavelength in SI units\n",
    "    print('wavelength: ' + str(wvlen));\n",
    "    ## =====================STRUCTURE======================##\n",
    "\n",
    "    ## Region I: reflected region (half space)\n",
    "    n1 = 1;#cmath.sqrt(-1)*1e-12; #apparently small complex perturbations are bad in Region 1, these shouldn't be necessary\n",
    "\n",
    "    ## Region 2; transmitted region\n",
    "    n2 = 1;\n",
    "\n",
    "    #from the kx_components given the indices and wvln\n",
    "    kx_array = k0*(n1*np.sin(theta) + indices*(lam0 / lattice_constant)); #0 is one of them, k0*lam0 = 2*pi\n",
    "    k_xi = kx_array;\n",
    "    ## IMPLEMENT SCALING: these are the fourier orders of the x-direction decomposition.\n",
    "    KX = np.diag(kx_array/k0);\n",
    "    KX2 = np.diag(np.power((k_xi/k0),2)); #singular since we have a n=0, m= 0 order and incidence is normal\n",
    "    #print(KX2)\n",
    "    KZ2 = ((n1)**2 - np.diag(KX2)).astype('complex');\n",
    "    KZ = np.sqrt(KZ2)\n",
    "#     print('KZ2: '+str(KZ2)+': imag Kz: '+str(KZ))\n",
    "#     KZ_mask = np.imag(KZ)<1e-2;\n",
    "#     KZ_mask = np.expand_dims(KZ_mask, axis = 1)\n",
    "#     print(KZ_mask.shape)\n",
    "    \n",
    "    ## construct matrix of Gamma^2 ('constant' term in ODE):\n",
    "    A = KX2 - E; #conditioning of this matrix is not bad, A SHOULD BE SYMMETRIC\n",
    "    #sum of a symmetric matrix and a diagonal matrix should be symmetric;\n",
    "    AO = np.block([[zeros, I],[A, zeros]])\n",
    "    beigenvals, bigW = LA.eig(AO);\n",
    "    print('conditioning of big block: '+str(np.linalg.cond(AO)))\n",
    "    ## SORTING IS REQUIRED\n",
    "    sq_beigenvals = np.square(beigenvals);\n",
    "    \n",
    "    #try rounding...\n",
    "    rounded_beigenvals = np.array([round(i,12) for i in beigenvals])\n",
    "    #print(rounded_beigenvals)\n",
    "    quadrant_sort = [1 if np.real(i)>=0 and np.imag(i)>=0 else 0 for i in rounded_beigenvals];\n",
    "    #print(quadrant_sort)\n",
    "    sorted_indices = np.nonzero(quadrant_sort)[0]\n",
    "    vsorted_indices = [i for i in range(len(quadrant_sort)) if i not in sorted_indices]\n",
    "    print(vsorted_indices)\n",
    "    print(sorted_indices)\n",
    "    #sorted_indices = np.argsort(np.real(rounded_beigenvals))\n",
    "    sorted_eigenmodes = bigW[:,sorted_indices];\n",
    "    Vsorted_eigenmodes = bigW[:, vsorted_indices];\n",
    "    #print(sorted_eigenmodes)\n",
    "    #adding real and imaginary parts seems to work...\n",
    "    sorted_eigenvals = beigenvals[sorted_indices]\n",
    "    Wp = sorted_eigenmodes[PQ:,0:PQ]\n",
    "    V = Vsorted_eigenmodes[0:PQ, 0:PQ ]\n",
    "    print(V.shape)\n",
    "\n",
    "    eigenvals_wp = (sorted_eigenvals[0:PQ]);\n",
    "\n",
    "    # when we calculate eigenvals, how do we know the eigenvals correspond to each particular fourier order?\n",
    "    eigenvals, W = LA.eig(A); #A should be symmetric or hermitian\n",
    "    #print('sorted eigenvals: '+str(sorted_eigenvals))\n",
    "\n",
    "#     plt.plot(sq_beigenvals);\n",
    "#     plt.plot(eigenvals)\n",
    "#     plt.show()\n",
    "    \n",
    "    W = Wp;\n",
    "    #we should be gauranteed that all eigenvals are REAL\n",
    "    eigenvals = eigenvals.astype('complex');\n",
    "    \n",
    "    #real parts match, but the imaginaries don't\n",
    "    Q = np.diag(eigenvals_wp);\n",
    "    #plt.plot(sorted(np.abs(beigenvals)),'.-');\n",
    "    colors = [(random.random(),random.random(),random.random()) for i in range(len(beigenvals))]\n",
    "#     plt.subplot(141)\n",
    "#     plt.plot(np.real(beigenvals), np.imag(beigenvals), '.'); plt.title('1st'); plt.ylim([-4,4])\n",
    "#     plt.subplot(142)\n",
    "#     plt.plot(np.real(eigenvals_wp), (np.imag(eigenvals_wp)), '.r', markersize = 10)\n",
    "#     plt.plot(np.real(np.sqrt(eigenvals)), abs(np.imag(np.sqrt(eigenvals))), '.g', markersize = 5)\n",
    "\n",
    "#     plt.ylim([-4,4])\n",
    "#     plt.subplot(143);\n",
    "#     plt.plot(np.real(eigenvals), np.imag(eigenvals), '.'); plt.title('2nd')\n",
    "#     plt.ylim([-4,4])\n",
    "#     plt.subplot(144)\n",
    "#     plt.plot(np.real(np.sqrt(eigenvals)), np.imag(np.sqrt(eigenvals)), '.g'); plt.title('sqrt 2nd')\n",
    "#     plt.ylim([-4,4])\n",
    "\n",
    "#     plt.show();\n",
    "#     plt.subplot(121)\n",
    "#     for b in beigenvals:\n",
    "#         plt.axhline(np.imag(b), color='r')\n",
    "#     for b in np.sqrt(eigenvals):\n",
    "#         plt.axhline(np.imag(b))\n",
    "\n",
    "#     plt.plot((np.imag(np.diag(Q))),'.g-', markersize = 10)\n",
    "#     plt.plot((np.imag(np.diag(Qn))), '.c-', markersize = 10)\n",
    "#     plt.title(str((len(beigenvals), len(np.diag(Q)))))\n",
    "    \n",
    "#     plt.subplot(122)\n",
    "    \n",
    "#     for b in beigenvals:\n",
    "#         plt.axhline(np.real(b), color = 'r')\n",
    "#     for b in np.sqrt(eigenvals):\n",
    "#         plt.axhline(np.real(b))\n",
    "        \n",
    "#     plt.plot((np.real(np.diag(Q))),'.g-', markersize = 10)\n",
    "#     plt.title('real')\n",
    "\n",
    "#     plt.plot((np.real(np.diag(Qn))), '.c-')\n",
    "#     plt.show()\n",
    "\n",
    "\n",
    "    V = W@Q; #H modes\n",
    "#     plt.subplot(121)\n",
    "#     plt.plot(np.real(Va))\n",
    "#     plt.subplot(122)\n",
    "#     plt.plot(np.real(V))\n",
    "#     plt.show();\n",
    "    \n",
    "    ## THIS ENFORCES a PARTICULAR SIGN...\n",
    "    X = np.diag(np.exp(-k0*np.diag(Q)*d)); #this is poorly conditioned because exponentiation\n",
    "    ## pointwise exponentiation vs exponentiating a matrix\n",
    "\n",
    "    ## observation: almost everything beyond this point is worse conditioned\n",
    "    k_I = k0**2*(n1**2 - (k_xi/k0)**2);                 #k_z in reflected region k_I,zi\n",
    "    k_II = k0**2*(n2**2 - (k_xi/k0)**2);   #k_z in transmitted region\n",
    "    k_I = k_I.astype('complex'); k_I = np.sqrt(k_I);\n",
    "    k_II = k_II.astype('complex'); k_II = np.sqrt(k_II);\n",
    "    Y_I = np.diag(k_I/k0);\n",
    "    Y_II = np.diag(k_II/k0);\n",
    "\n",
    "\n",
    "    delta_i0 = np.zeros((len(kx_array),1));\n",
    "    delta_i0[num_ord] = 1;\n",
    "    n_delta_i0 = delta_i0*j*n1*np.cos(theta); #this is a VECTOR\n",
    "\n",
    "    ## design auxiliary variables: SEE derivation in notebooks: RCWA_note.ipynb\n",
    "    # we want to design the computation to avoid operating with X, particularly with inverses\n",
    "    # since X is the worst conditioned thing\n",
    "\n",
    "    Wi = np.linalg.inv(W);\n",
    "    Vi = np.linalg.inv(V);\n",
    "    Oi = 0.5*np.block([[Wi, Vi],[Wi, -Vi]])\n",
    "    f = I;\n",
    "    g = j*Y_II; #all matrices\n",
    "    fg = np.concatenate((f,g),axis = 0)\n",
    "    #ab = np.matmul(np.linalg.inv(O),fg);\n",
    "    # ab = np.matmul(Oi, fg);\n",
    "    # a = ab[0:PQ,:];\n",
    "    # b = ab[PQ:,:];\n",
    "\n",
    "    a = 0.5*(Wi+j*Vi@Y_II);\n",
    "    b = 0.5*(Wi-j*Vi@Y_II);\n",
    "    fbiX = np.matmul(np.linalg.inv(b),X)\n",
    "\n",
    "    #altTerm = (a@X@X@b); #not well conditioned and I-altTermis is also poorly conditioned.\n",
    "    #print(np.linalg.cond(I-np.linalg.inv(altTerm)))\n",
    "    #print(np.linalg.cond(X@b)); #not well conditioned.\n",
    "\n",
    "    term = X@a@fbiX; # THIS IS SHITTILY CONDITIONED\n",
    "    # print((np.linalg.cond(X), np.linalg.cond(term)))\n",
    "    # print(np.linalg.cond(I+term)); #but this is EXTREMELY WELL CONDITIONED.\n",
    "    f = np.matmul(W, I+term);\n",
    "    g = np.matmul(V,-I+term);\n",
    "    T = np.linalg.inv(j*np.matmul(Y_I,f)+g);\n",
    "    T = np.matmul(T,(np.matmul(j*Y_I,delta_i0)+n_delta_i0));\n",
    "    R = np.matmul(f,T)-delta_i0;\n",
    "    T = np.matmul(fbiX, T)\n",
    "\n",
    "    ## calculate diffraction efficiencies\n",
    "    #I would expect this number to be real...\n",
    "    DE_ri = R*np.conj(R)*np.real(np.expand_dims(k_I,1))/(k0*n1*np.cos(theta));\n",
    "    DE_ti = T*np.conj(T)*np.real(np.expand_dims(k_II,1))/(k0*n1*np.cos(theta));\n",
    "    #print(np.sum(DE_ri))\n",
    "\n",
    "    #print(np.sum(DE_ri))\n",
    "    spectra.append(np.sum(DE_ri)); #spectra_T.append(T);\n",
    "    spectra_T.append(np.sum(DE_ti))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd8m9X1/99X3jueWU4cj2xnQJxF\nIIMRMiEN0BJmQlPotwtaoNDBKC2l9EcXpS1lBihNCAUCBMgCQsgeJIQszzixEzu27Njxnvf3x5Uc\nxZFsydaw5Pt+vfySpWfcY1l6znPPOfdzhJQSjUaj0WgADJ42QKPRaDQ9B+0UNBqNRtOGdgoajUaj\naUM7BY1Go9G0oZ2CRqPRaNrQTkGj0Wg0bWinoEEIsUII8TtP26GxjRDid0IIoxCi2NO2eAtCiGoh\nRIqn7fA2tFPQdAshhBRC1Ji+gKeEEH8WQvi12+daIcQWIUSVEKJUCPGFEOI6IYS/6bhJFvveajpn\n+9eOtTvnUtN+33b9X2k/QojHhRD/cfI5BwH3A6OklP2ceW47x58lhPhcCFEphMi3st3yM1AmhPhU\nCPEdN9u4WQix3PI1KWW4lDLPnXb4AtopaJzBOCllODAD+A5wl3mDEOJG4G3gdSAR6As8CiyUUjYD\nO0zHmZkOHLPy2pZ2Y94JlJsevQahcPR7lwSUSSlLujCev6PHWKEGeAV4sIN9zJ+B4cAK4DkhxGNO\nGNtZf4PGXqSU+seDP8BDwCmgCsgErjK9Pgl1wawAioDngECL4yTwAyDbdOxvgVTTMeeA1eb9gZlA\nIfBLwAjkA7danGsF8DuL5wuAA6axtwNjO7BfAmkWz1cD/zD9LoCTwIMdHP8I8KHF8yPAUiuv3Wbx\nPAloBW4AmoG+HZw/DfgCqDT97W+1s/0nQJ5p2/8DDBbb7wKOAmeB9UCSxbbRwEaUYzpjem/nAI1A\nE1ANfG3adzPwJLANqDPZtMx07irT+PfYsP9q0zGtpnOuML1+HXDY9D/aDIy0OCbf9Lk6CDQA/lbO\n+zegwPRZ2QdcYcdn9Wogv7PPgOm1G4F6INbGuS4F9pv+/reBt8yfQc5/Xh8CioE3gGhgLVBq+n+s\nBRJN+z8JtJjGqwaea28X6jP+D+Aj05i7gFQLe2ajvn+VwD9Nn5nlnr4+eOLH4wb05h/UXVUBMMD0\nfIj5gwpMAKYA/qbXjwL3WRwrgQ+ASNMFqgH4FEgBolAX0jtN+85EXTz/DASh7sJrgOGm7SssvpCX\nAiXAZMAPdSeeDwTZ+Bssv3gjUA7spxbPJZDcwXswA3VhNQBxwAkgFHWhNb/WCgy2OOYRYLfp92+A\nn3Vw/pXAr0znCgYub2f750AMMBjIMl8IgEVADjDS9D/4NbDdtC3C9HfebzpnBDDZtO1x4D/tbNiM\nco6jTecKAOajnLgwvQe1wKU2/oaZQKHF82Gm/981pnP93GSr+SYgH+XUBwEhNs55GxBrsud+1MU3\nuJPPqyNOIcD0mZtrZf9A0//5XtN+i1HO1NIpNANPoz6vISZbbzB9NiJQjmRNu/d4uS27UJ/xctTN\nlj/wJrDKtC0O5RwXm7bdi3Ls2inoHze/+eqOscT0ZQvoZN/7gPcsnktgmsXzfcBDFs//BPzV9Lv5\nSxZmsX018Ijp9xUWX8h/Ab9tN3YmMMOGXdL0haox/b4SkwMBpples3mxQV1U64FxwLeAN02v77R4\n7Xi7Y7IxOUjgF5juyG2c/3XgBUx3lVZsn2Px/AfAp6bfPwG+a7HNgLpwJwFLgP02xnsc607hiU7+\nv2uAe21sm8mFTuERYHU7204BM03P84G7HPwsnkWFgDrax26nYHq9GIsZqcXr0032CovXtnKhU2js\n5HMzHjjb7j3uzCm8ZLFtHnDM9PsdwA6LbQJ1s9YrnYLOKXgQKWUO6mL/OFAihFglhBgAIIQYJoRY\nK4QoFkKcA36PuqOx5IzF73VWnodbPD8rpayxeH4CGGDFrCTgfiFEhfkHdcdpbV8zl5rG+g5qhhFm\ner3M9Njf1oFSynpgN+pCMR340rRpq8VrbfkEIcQ0IBlYZXrpv8AYIcR4G0P8HPUl3y2EOCyEuKvd\n9gKL3y3fkyTgbxbvQbnpPANR70eurb/JBpbjIISYK4TYKYQoN51/Hhf/f20xwGQrAFLKVtP5B9oa\nrz1CiPuFEEdNyeMK1OzS3vE7RQgRAMSj3rf2DABOSdMV2Ia9pabPhvl8oUKIfwshTpi+D1uAPu2L\nGjrBsnKrlvPfjwGW45vsKnTgvD6FdgoeRkr5Xynl5aiLkERNmUHdsR8DhkopI1Exa9GNoaKFEGEW\nzwcDp63sVwA8KaXsY/ETKqVc2cnfIaWUq1E5jUdNL2eazndDJ7ZtQV38r+C8U/jS4jXLJPOdqPfh\ngKk8c5fp9Tts2FUspfyelHIAcA/wTyFEmsUugyx+t3xPClBxfsv3IURKud20LdXG3yI7e10IEQS8\nAzyDyof0AT7G/v/vadTnxXw+Yfo7TtlhB0KIK1Dx+m8D0abxKx0Y3x6uR81Od1vZVgQMNNltZlC7\nfdrbfz8q3DrZ9H2Ybnpd2NjfEYpQRRDqhMquRNu7+zbaKXgQIcRwIcSVpotEPeruvsW0OQIVlqkW\nQowA/s8JQ/5GCBFouigsQMVl2/Mi8H0hxGRTpUyYEGK+ECLCzjH+ANwthOhnuuP6GfCIEGKZECJS\nCGEQQlwuhHjB4pgtwCzUheGI6bWtqDDCeNN2hBDBqAvZ3abXzT8/Bm61VqUihLhJCGH+gp9FXTxa\nLHZ5UAgRbSr7vBeV8AR4HviFEGK06TxRQoibTNvWAv2EEPcJIYKEEBFCiMmmbWeAIZ1UGAWiYuWl\nQLMQYi4q0Wkvq4H5QoirTHfk96NyStvtPD4CdcEuBfyFEI+iclNWMf3PglHxfyGECBZCBNrYN0YI\ncSsqqfu0lLLMym47UP+DH5nKkq9Hxfo7s7kOqBBCxACPtdt+BpVP6wofoWabi0yfoR8Cbi/97Slo\np+BZglAXUSNqapuAmhEAPADcgqqUeJHzF6uuUoy6KJ5GJdm+L6U81n4nKeVe4HuoaqezqATmUnsH\nkVJ+g6rceND0/H+cL1M9jfry/g543+Kw7ajwxS5zSMF0MSkFSqSU2ab9FqEuDK+bZgDFUspi4GVU\nUnyOFZMmAruEENWoxPy9UsrjFtvfR+VjDqAuDi+bxn8PNWtbZQpXHALmmrZVoZK8C1HvazbKqcF5\nR1smhPjKxntUhap6Wo16j28x2WYXUspMVKL476jPzkJUiW+jnadYj8qZZKHCUPV0HG6ajnrfP0bN\npuqADe32+dr0HucAy1HFBo9iBZOdi4HvoqqnbkM52oYObPgrKuFsROWb1rXb/jfgRiHEWSHEsx2c\nx5o9RuAm4I+okOcoYG8n9vgs4sKwnsYXEULMRCU/e+2U2BpCCIkKz+V42pbejhBiF/C8lPLVHmCL\nAZVTuFVK+bmn7XE3eqag0WjcjhBihhCinyl8dCcwlovv/t1pz7VCiD6mUK45f7fTU/Z4Epc5BSHE\nK0KIEiHEIRvbhRDiWSFEjhDioBDiUlfZotFoehzDga9RCe77gRullEUetGcqqqLMHI5bJKWs86A9\nHsNl4SMhxHTU6sLXpZTpVrbPQyUI56HKGP8mpZzcfj+NRqPRuA+XzRSklFuwXqNs5nqUw5BSyp2o\nmmOb9ewajUajcT2eFJoayIUVD4Wm1y6aQgoh7kaVIRIWFjZhxIgRbjFQ40M01YIxB+KHg3+Qp63R\neBuNNVBZCAY/CAg1/YR41Wdp3759RillfGf7edIpWFsoYzWWJaV8ASVVQEZGhty7d68r7dL4Ituf\ngw2/givvhOkPeNoajbcgJex5CdY9DJFDILgPlByFVpM4QFAU9B8LA8ZDf9NPTAoYel4NjxDiROd7\nedYpFHLhKsZErK+w1Wi6jzFTPWZ+rJ2Cxj6a6uGjn8GBN2HotbD4BQjpA80NUHIEir6G0wfU464X\noMW0rCEwQjmK/uOh/zjlMGLT1CzDC/CkU/gAtaJxFSrRXOnh6gONL1OapR5P7YNzRRCp01eaDqgo\ngNW3w+n9MOMhmPHw+bt//yAYcIn6mWDav6VJzSCKvoYik6PY+zI0m+SbAsKUg5jxc0idZXXInoLL\nnIIQYiVKpiBOCFGIWpYeACClfB61OnIeagVkLUpfXqNxPlJC6TEYPBVO7oCsdZChP24aGxzfAm8v\nheZGuHkljJjX+TF+AabZwVjgdvVaSzMYs5STOH0AsjfAG9+CK38Fl9/fI0NM4EKnIKVc0sl2idIY\n0WhcS00p1FfAqOuhqkiFkLRT0LRHStjxD9j4qAr33PwmxA3t+vn8/KHvKPUz/hZoeBQ+/Al89jso\n3Aff+heERDvPfifRM12VRuNMSk35hLhhMHwe5H0BDdWetUnTs2ishXeWq2KEEfPge592zyFYIygc\nbngZ5v4RcjbCCzOh6KBzx3AC2ilofB9zkjl+hHIKLQ2Q+5lnbdL0HMqPw8vXwKF34KpH4dtvQJC9\nosAOIgRMvgeWfaLCUy9fA/vfdM1YXUQ7BY3vU5qpKkIiB6i8QnAfyPzE01ZpegI5m9Qde2UB3Po/\nuOJ+deF2NYMmwT1bIHEivP8D+OAnqtqpB6Cdgsb3Kc1UoQAhVJx32LUq2dzS7GnLNJ5CSvjyT/Cf\nGyFyINy9GYZe7V4bwuPh9jVw+U/hq9fglWvhrF1LCVyKdgoa38eYpVYymxk+D+rKodBaUzCNz9NQ\npcpNP30C0hfD8o1qwZkn8POHqx+Hm/8L5XnwwgzI3uQZW0xop6DxbeorVcWRpVNIuwr8AuHYR56z\nS+MZjNnw4lXqfz/7SZX4DQzr/DhXM2K+mq1EDoQ3b4TNf4DWVo+Yop2Cxrcxmpq2xVk4haAISJ6u\nSlN1k6neQ+Yn8OKVUGtUYZvLfuSe/IG9xKbCdzfC2O/A5qfgvzdBbUeaoq5BOwWNW9iWY2TCbzdS\nWuXmDoelpo6jljMFgOFz1XTdmOVeezSeoWA3rFyiwkR3fwEpMzxtkXUCQ+Fbz8P8P6tFdP+eoVZV\nuxHtFDRu4aUv8yiraeTQ6Ur3DlyaqUJFfZIufH3YXPWoQ0i+T2srfPIQhPeFpWuhz6DOj/EkQsDE\n78KydSBb4eXZsG+F22a12iloXM7pijq+yCoFIK+0xr2DG7PU6lS/dov3owYqwTJdmur7HHwLTn+l\nErquWn/gChInqLLVpGnw4b3w/o+gyfXN4LRT6EW0tHomfv723kJaJQQHGMgtdfNK4tLMi0NHZkbM\nh8I9UF3iXps07qOhGjY9DgMnqFi9txEWC7e9A9N/Dgf+Azv/5fIhtVPoJXyeWcL432ygoLzWreO2\ntEpW7y3giqFxjOofSZ47nUJTHVScuDDJbMnwuYBUaxY0vsnWv0B1Mcz5Q48VoOsUg58S0Vv6EUz9\nkeuHc/kImh7B3vxyqhqaeX1HvlvH3Zpj5FRFHd+ZOIjU+HD3ho/KclRMNn6Y9e190yFqMBz72H02\nadzH2ROw/e8w5ia1gtjbGXI5+Ae6fBjtFHoJ5ovxqj0F1DS4byXvW3tOEhMWyDWj+pISH05JVQNV\n9U3uGbxNCM/GTEEIJX6W97lqt6jxLTY+qu6yr/6Npy3xKrRT6CXkllYzsE8IVfXNvPNVoVvGNFY3\nsPHIGRZfMpAgfz9S4tUiIbfNFoxZIAwq0WyL4XNVI5S8ze6xSeMe8rfBkTUw7T5VVKCxG+0UegEt\nrZJ8Yy0LxvZnXGIUK7bl0+qGpPM7+wppapHcPEmVAKbGhwOQZ3RTXqE0E6KHQECw7X2Spqk+uzqE\n5Du0tsC6hyAyES77saet8Tq0U+gFFJ6tpbGlldT4cJZNSybPWMMX2aUuHVNKyVt7CshIiiYtQZUB\nDo4Jxc8gyC1x40zBVujIjF8ADL1GJZtbW9xjl8a17P8PFH8D1/xGLQbTOIR2Cr0AcxloakIY88b0\nJyEiiFe2HnfpmLuPl5NnrOHmSYPbXgv0NzA4JtQ9M4WWZiVxYSvJbMmIeUr6oHCP6+3SuJb6c/DZ\nb2HQFEi/wdPWeCXaKfQCzDH8lLhwAv0N3D4liS+zjeSUVLlszLf2FBAR5M+8Mf0ueD01Psw9OYWz\n+dDa1PlMASDtajAEKC0kjXez5f+p9qtznupZukZehHYKvYDc0mpiwgKJDlPlbLdMHkygv4FXt+W7\nZLzK2iY++qaI6y8ZQGjghSuJU+LDyTPWuH4hnWW3tc4IjlLlfjqv4N2U5arFXeNvhYGXetoar0U7\nhV5AbmkNKXHn5YFjw4NYNH4A73xVSEVto9PHe//rUzQ0t3LzxMEXbUuJC6OxuZXTFS5ert9Wjmpn\nn90R86Es+7yqqsb72PAI+AeplpqaLqOdQi8gr7S6rfLHzLJpydQ3tbJqT4FTx5JSsnJ3AekDI0kf\nGHXR9tQEZYfL5S5KMyFiAARH2rf/sDnqUYeQvJO8zZD5EVzxM4jo1+nuGttop+DjVNY2YaxubFsj\nYGZk/0impMTw+vZ8mluc18zjm1OVHC06Z3WWALTNWHJdnVcwZtqXZDbTZxD0G6tDSN5ISzOs+4VS\nwp3yQ09b4/Vop+Dj5JoqfdrPFADumpbM6cp61h8+47TxVu4uICTAj+vGD7C6PSYskKiQANdqIEmp\nwkD2JJktGT4PCnZBjdE1dmlcw1croOQIzP5dx2tSNHahnYKP01Z5FH9xy8GrRvZlUEwIr25zTnlq\nTUMzHxw4xfyx/YkMDrC6jxCC1Pgw14aPzp2Cxmrb6qi2GDEPLZDnZdSdhc+ehCFXwMiFnrbGJ9BO\nwcfJLa0mwE8wKObiRTx+BsGdU4ew98RZDhZWdHusjw4WUdPYws0TO25ikuJqYTxzktlRp9BvrFoF\nq3sseA9f/BHqK3QJqhPRTsHHySutZnBMKAF+1v/V3544iLBAP6eUp67ac5K0hHAmJEV3uF9KfJhr\nhfE6E8KzhRBKCyn3M7c0M9F0k9Is2P0CXHoH9BvjaWt8Bu0UfJzc0hqr+QQzkcEB3JQxiLUHT1Ny\nrr7L42SdqeKrkxXcPHEQopM7tjYNJFfNFoyZEBINYXGOHztiHjTVQt4XzrdL41zW/xICQmHWrz1t\niU+hnYK7yP0cXrzSrRLNzS2tnCiraSsDtcWdlw2huVXyn50nujzWqt0FBPgJFl+a2Om+qWa1VFfJ\nXZSaNI+6Ek5IuhyCIlV5o6bnkr0RcjbCjJ9DeLynrfEpeo9TyP0M1vzAbc2vL2LPS3BqH5zY4bYh\nC87W0dQiL1i4Zo3kuDCuHJ7Am7tOUt/kuChcfVML7+4vZPbofsSEdd4EZHBMGH4G4dqZgqP5BDP+\ngUr2InOdaviu6Xm0NKlZQkwqTLrH09b4HL3HKVSeggNvqgbe7qahSt3ZABx3X1git8QshNfxTAHU\nYraymkY+/Pq0w+OsP1xMRW0TS2ysTWiPWRjPJRVINUaoLeu6UwBVmlpTopy4puex5yWlgHvtk27p\nRNbb6D1OYeQCJXp26F33j525DloaICQGjm9x27Dm8ExqXOdOYVpaLMP6hvPKtnykg7Opt/YUMCgm\nhMtSY+0+JiXORcJ4XU0yWzL0ajD46xBST6SmDDY/BSmzzq9C1ziV3uMUQqIh7So4vMb9YYEjayCi\nP0z6HhR9rWqr3UBuSQ1x4YFEhVpfM2CJEIJl05I5WnSOXcfL7R7jRFkN23PL+E7GIAwG+2P4KfFh\nHHeFMF6bEJ4Dq5nbExINSZfp0tSeyObfQ0O1LkF1Ib3HKQCMXgznCqFwt/vGNIeORl6n7m6QkL/V\nLUPnGatJsWOWYGbR+IH0CQ1waDHbW3sKMAi4KaPjtQntSY0Pp8EVwnilWRAQptYbdIfh86H0mFLe\n1PQMzhyBva/AxO9CwkhPW+Oz9C6nMHwu+AW5N4SUtV6FjkYvgoETVAmdm0JIuaU1pCZ0nGS2JCTQ\njyWTBrPhyBkKyms73b+ppZW39xVy5YgE+kY6Ji+QEu8iYTxjplJGNXTzoz18rnrUs4WegZSw/heq\nMmzmLzxtjU/jUqcghJgjhMgUQuQIIR62sn2wEOJzIcR+IcRBIcQ8V9pDcCQMm63COe5qvXj4PQjv\npzpB+QfC4KlucQpnaxopr2l0aKYAcMfUJAxC8Nr2/E73/fxYCaVVDTbF7zqirSzV2XmF0qzuJZnN\nRCdB33StmtpTyPtcKaHO/AWExnjaGp/GZU5BCOEH/AOYC4wClgghRrXb7dfAainlJcDNwD9dZU8b\noxdD9Rk4sc3lQ9FQDTmbYNT15+9ck6ersESV80TorNGWZHZgpgDQPyqEuen9eGtvAdUNzR3uu2pP\nAX0jg5g53PE6cbMwnlNnCg1VKjwY1418giXD58HJHVBrf45F4yK2PAORAyFjmact8XlcOVOYBORI\nKfOklI3AKuD6dvtIwCx4HwU4Xg/pKMOuVTHnQ++4fCiy1kFzvQodmUmerh5dPFvItWjB6SjLpiVT\nVd/MO/sKbe5TVFnH5swSbpowCH8bEhodIYQgxdmtOY1Z6tEZMwVQISTZqkKAGs+Rv03dxE27VzXR\n0bgUVzqFgYBlB5dC02uWPA7cJoQoBD4GfmztREKIu4UQe4UQe0tLS7tnVWAYDJ8DRz5Qi2BcyZE1\n50NHZvqPU+0fXbxeIbe0mkA/A4nRIQ4fe+ngPowb1IcV2/NptVEd9PbeQlolfNvBBLMlqfHhzl3V\nXGp2Cna04LSHAZeoqjEdQvIsW/4fhMUrjSONy3GlU7BWL9b+CrMEWCGlTATmAW8IIS6ySUr5gpQy\nQ0qZER/vhCXt6TdAXblrL8wN1arqaNR1FyY9DX5K5tfVM4WSGpJiQ7t8F3/XtCEcN9bwRdbFTri1\nVfLWngIuT4tjcOzF6qv2khIfxplzThTGM2aqtSjRyc45n1kgL+dTaOq6LpSmGxTuVfmEy34MAY7f\n4Ggcx5VOoRCwvI1M5OLw0HeB1QBSyh1AMNAFFTMHSbtaVTG4sgope70KHY1adPG25OlQcQLO5rts\n+DzjxS04HWFuen8SIoJ4xUp56tYcI6cq6rh5UtdnCXA+tHXc6KQQUmkmxKaCn79zzgeqNLWpxq2L\nDjUWbHlGrRvJuMvTlvQaXOkU9gBDhRDJQohAVCL5g3b7nASuAhBCjEQ5hW7Gh+zAPwhGLICja6G5\nwTVjHF4D4X1h8JSLt7XlFb50ydBNLa2cLKu12ljHXgL9DdwxNYkvs41kn6m6YNtbewqIDg3gmlF9\nu2VnWoKTK5BKM52XZDaTfAUEhusQkicoOghZn6gWm0ERnram1+AypyClbAZ+BKwHjqKqjA4LIZ4Q\nQlxn2u1+4HtCiK+BlcBS6ajGQldJXwwNlSo04Gwaa84vWDP4Xbw9fgSEJbgsfHWyvJbmVtmtmQLA\nkkmDCfQ38KpFeWpZdQMbjhRzw6WJBPlb+dscwCyM55QKpOYGOHvceUlmM/5BaiV85idaIM/dfPkn\nNaOf9D1PW9KrcOk6BSnlx1LKYVLKVCnlk6bXHpVSfmD6/YiUcpqUcpyUcryUcoMr7bmAlJlqWnrY\nBSGkrPXQXHdh1ZElQqjZwvEtLlFtNQvhdWemABAbHsSi8QN496tCKmobAXj3q1M0tchuh45AzUYG\nRYc4Z6ZQlqsqhZyVZLZk+HyoLoai/c4/t8Y6pZlw5H2YdDeE9PG0Nb2K3rWi2RK/AHUnf+xjaOx8\n9a5DHFmjZgKDp9reJ3m6Wi9hLqN0InlGc1/m7s0UQJWn1je1snJ3AVJKVu45SUZSNGkJzpnOp8aH\nO2emYNY8cnb4CGDoNSD81GdF4x6+/JNKLE/5gact6XX0XqcAKoTUVAPZTpygNNZA1gZT1VEH4RUX\nrlfILakmLjyIqJDOhfA6Y2T/SKamxPL6jnx25JWRV1rDdzrpwewIZmE8W6WvdlOaCQglceFsQmO0\nQJ47Kc+Db95WyeUw+5V3Nc6hdzuFIVeoO3pnhpCyN6jQkbWqI0uih0DUYLV038nkGWvaZCScwbJp\nQyiqrOfBtw8SEeTP/LH9nXbuFJMw3qnuCuOVZkKfwa4rWxw+F0oOQ7n9YoGaLrL1L6q0+LKfeNqS\nXknvdgoGPyVBkbVeSSQ4g8Om0FHSZR3vJwSkTFeKqU7WYcotrbarsY69XDWyL4NiQjhVUcd14wcQ\nGui8ks+2fs3dLUs1Zrkmn2BmxHz1eEz3WHApFQVwYCVMuBMiulfdpukavdspgFrI1lyvGuF0l8Ya\n5WBGLuw4dGQmeQbUV0DxN90f20R5TSMVtU2dtuB0BD+D4K5pyQihKpKciTkZbk6Od4nWFjBmd6+H\nQmdED4G+Y+DYWteNoYFtf1OPepbgMbRTGDQZIgY4RwvJHDqyVXXUniFXqEcn5hXMSVtnzhQA7pw6\nhE0/m0H6wCinnjc2LJDIYP/uyV1UnFDy5N3ptmYPIxfAyZ1QXeLacXorVcXw1eswfgn0cV7eSuMY\n2ikYDCrhnLOp+x3RDq9RGi1J0+zbP7K/qpZxolPIK7W/BacjGAyi2+serCGEIDUhvHtlqeYWnM5e\no9CekQsBqUNIrmL736G1CS7/qact6dVopwBKTru1qXtf9sZaNVOwN3RkJnk6nNgOzY1dH9uC3NIa\nAv0NDOyCEJ6nSInrZllqqQvLUS1JGKV0lXQIyfnUlKmuamNugpgUT1vTq9FOAWDgpdAnqXtaSNkb\noKm286qj9iTPUGWxp7/q+tgW5JZUkxyrVgp7C2ZhvM76N9jEmKXUaF29yEkIFULK+wLqK107Vm9j\n5z+hqQ6uuN/TlvR6tFMA9WVPX6zKQ2vKunaOI2sgNM7+0JGZIZcDwmkhpDyjYy04ewJtFUhdnS2U\nZro2yWzJyOvUrDLLfYvvfZ66Ctj9gqoEdHUIUNMp2imYGb0YZAscfd/xYxtrz1cdOarQGRoD/cY4\nxSk0Nrdysry2S411PEm3WnNKaRLCc9PFZGCGmpUc+9A94/UGdr8ADedg+gOetkSDdgrn6TcGYod2\nLYSUs1GFjuytOmpP8nQo2NVtuY2T5TW0tEqvmykMjg3FILo4U6gqgsYq991hGgwwYh5kb1LhDk33\naKhSoaNhc9V3UONxtFMwYw7tw8LSAAAgAElEQVQh5W9VpXGOcHgNhMZC0uVdGztlJrQ0KsfQDXJK\nut6C05ME+fsxOCa0rYWoQ7ir8siSkQtVHij3c/eN6avsfUVV/elZQo9BOwVLRi8GpFJntJemuq6H\njswMngIG/26HkMy1/t1VR/UEKV0VxjMLCrorfARqfUlwlK5C6i5NdaoMNWUWJGZ42hqNCe0ULEkY\nAQmjHVvIlr1R3TU6WnVkSVAEDJzQbaeQW1JDQkQQEcHdF8JzNylxXRTGK81UF+jwBNcYZg2/ABg2\nRzXeaelixZRGLVSrKYUZP/e0JRoLtFNoT/piFcapLLRv/yOm0JF5dXJXSZ6uylK7UerY3RacniQ1\noYvCeOYks3BzCe6IBSrscWKbe8f1FZobYOtfVbVeZzphGreinUJ70herx8Pvdb5vU53STOpO6MhM\n8nTVJObE9i4dLqUkt6TaK0NHQJtWk8PCeEY3lqNaknYV+IfoEFJXOfBfqDqtcwk9EO0U2hOTAgMu\nsS+ElLOp+6EjM4mTwD+4yyGksppGztU3e+1MIaUraxVqy1X4wZXqqLYIDFOO4eha3abTUVqaYOuf\nVcg0ZZanrdG0QzsFa4xeDKf3q2YfHXHYSaEjgIBgJc7XRafgrBacniIuXAnjOZRs9kSS2ZIRC9Td\n7mndptMhvvkfVJyE6Q+6P+yn6RTtFKwx+lvqsaM1C011kLVOXRi6GzoykzwdzhyCGqPDh5rDLt46\nUxBCkBLvoDBeWzmqB8JHAMOuVVVjeiGb/bS2qFabfceoZL2mx6GdgjX6DFJ37R3lFXI+hcbqri9Y\ns0byDPWY/6XDh+aWVBPkb2BgH+8RwmtPSnyY407BP0R1sPMEoTFKpuSozivYzZH3oSxb5RL0LKFH\nop2CLUYvVnft5rvR9hxZAyExMGS688YccAkERijBNQfJLa0mOS4MgxcJ4bUnNT6c4nP19gvjGTMh\nLk2tMvYUIxaoi5ytz4nmPK2tsOUZFe4beZ2nrdHYQDsFW4xeBAjrIaSmOtXEfaQTQ0egzjVkWpfy\nCkoIzztDR2bMGkjH7Z0tlLq4Bac9jFigHo9+4Fk7vIGsT1Sf6yvu96wj13SI/s/YIqKfCg0cfleJ\nrlliDh05o+qoPcnToTzX/nUSQENzCwXltaQ6sQWnJ2irQLKnC1tjDVSe9FyS2Uxkf0icqENInSEl\nbPl/qq1p+g2etkbTAdopdET6YlXhcubQha8fWQMh0eoC7mzM5zxuf17hRFktrdL5LTjdTZJJGM+u\nfs3GbPXoqSSzJSMWQNEB1XReY51ja1WV1uU/c+7sWuN0tFPoiJHXg/C7MITUVK8WrI1YoOQOnE3C\naJWrOG5/XqGtHNXLhPDaE+Tvx6CYUHLtWcDW1m2tB+jvj1yoHvVCNus01cH6X0L8SBh/q6et0XSC\ndgodERarFEwPvXM+hJT7qZJqdmbVkSUGAyRfofIK7cNWNjCXo3rrGgVLUu0tSzVmqnLQntC6MTZV\nterUISTrbH9OrUuY+7SeJXgB2il0RvpiqDhxvl3mYXPoaIbrxkyeAedOdb54zkRuSTX9IoMJC/L+\nL5wSxqvuXBivNFM5BP9A9xjWGSMWwMntXVpj4tNUFqp1CaOuhxQXfmc0TkM7hc4YMR8MASqE1FSv\nqo5GzHdN6MiM2eHYGULK9cIWnLZIiQ+nvqmV05WdCOMZsyCuB+QTzIxcoLSrMj/2tCU9i42PAhJm\n/87TlmjsRDuFzgiJhrSr1UI2c+ho1LdcO2ZsKkQMsKs0VUpJXkm11+cTzNjVmrO5Ecpye1Y/335j\noc9gHUKyJH+bCr1Ou0+9NxqvQDsFe0hfrMI5m34DwX1cPw0WQlUhHd/SqdhaaXUDVQ3NbRdTb8dc\nltqhBlJ5nuqn3ROSzGaEgBELIe9z1WKyt9PSDJ/8HKIGwbR7PW2NxgG0U7CH4XOVgqkx03VVR+1J\nmQG1ZVBypMPdcs0tOL1U86g9ceGBRAT7dzxTMHqgBac9jFyg2qpmb/C0JZ7nqxWqlHv27yAw1NPW\naBxAOwV7CIqAobPV766qOmqPWXm1kxCSeaGXt69RMCOEUBVIHS1gKzWrow51j1H2MmgyhMbpEFJt\nOXz2O/UZHnW9p63ROIh2CvZy2U+UHpIrq44s6TNIVdd04hRyS2oIDjDQPzLYPXa5gZT4sLYZkFWM\nmUoEL7CHhcwMfqoIIXuDKkrorXz+e9VBcO7TWvTOC3GpUxBCzBFCZAohcoQQD9vY59tCiCNCiMNC\niP+60p5uMWgi3PSqe0sgk6erdo8d9AHOLVVJZm8WwmuPWRivxpYwXumxnrGS2RojFyoJFAcWH/oU\nxYdg78swcTn0He1pazRdwGVOQQjhB/wDmAuMApYIIUa122co8AtgmpRyNHCfq+zxSpKnQ8M5JaFg\ngzyj97bgtEWbMJ61lc2trWDM6VlJZkuSpyul26O9sMeClPDJQ6oYY+YvPG2Npou4cqYwCciRUuZJ\nKRuBVUD7AOP3gH9IKc8CSClLXGiP92GW5bZx11nf1ELh2Tqvbaxjiw4rkCpPQnNdz0sym/EPUs13\nMj9WDWV6E4ffgxNb4cpfq14TGq/ElU5hIGCpEFZoes2SYcAwIcQ2IcROIYTVVkxCiLuFEHuFEHtL\nS0tdZG4PJDxeaSHZyCvkl9UgpW/IW1jSJoxnrQLJnGTuqU4BVBVSbRmc3OFpS9xHYy1seER1VJuw\n1NPWaLqBK52CtSB3e+0Cf2AoMBNYArwkhOhz0UFSviClzJBSZsTHxzvd0B5N8nQ4uROaGy7aZE7G\n+tpMwSyMl2dtpmAuR+1Jq5nbk3YN+AX1riqkbX+Fc4Uw748q4a7xWlzpFAqBQRbPE4HTVvZ5X0rZ\nJKU8DmSinITGTPJ0aK6Hgt0XbTJfNH1tpgBKA8n6TOEYhMX37PBEUDikXqlUU+0UNfRqzubD1r9C\n+o2QdJmnrdF0E1c6hT3AUCFEshAiELgZaN+eag0wC0AIEYcKJ9mnAtdbSLoMhMFqCCm3tJoBUcGE\nBnq/EF57UuLDrQvj9YRua/YwcgFUFkDR1562xPVs+LWaHVzzhKct0TgBlzkFKWUz8CNgPXAUWC2l\nPCyEeEIIYW7Quh4oE0IcAT4HHpRSlrnKJq8kpI/q3WzFKfhCC05bpFoTxpPS1Je5B4eOzAybq5y5\nr1ch5W1Wf+MVP4Oo9ilDjTfi0nUKUsqPpZTDpJSpUsonTa89KqX8wPS7lFL+TEo5Sko5Rkq5ypX2\neC3J0+HUXmg4H2OXUpJbUk2Kl7fgtEWKNWG86hK1KKonJ5nNhMVC0jTfbrzT0gSfPAx9kmDqjz1t\njcZJ+F7cwRdJng5b/6KqWYZeA0BJVQM1jS0+O1M47xSqmT7MVFzgDUlmS0YuVKJwxuyeJ8nhDPa8\nDKVH4eb/0oQfhcePU1/fi1dy9xCCg4NJTEwkIKBrGm3aKXgDg6aAX6Bar2ByCr7SgtMW8eFBRAT7\nX5hsLu2hQni2GDFfOQVzeMWXqDEqOYvUK2H4PArz84mIiGDIkCEILW3hMaSUlJWVUVhYSHJycpfO\nobWPvIHAUEicdEFewdzH2Fea67RHCEFKe2G80kwIioSI/p4zzBGiElU+yBdDSJ/9FppqYM4fQAjq\n6+uJjY3VDsHDCCGIjY3t1oxNOwVvIXk6FB2EugpAzRRCA/3o50NCeO1JjQu7MKdgTjJ704Vn5EI4\ntQ8qT3X7VPVNLcieUOJ6+gDsew0m3XPBrE07hJ5Bd/8P2il4C/3HAbIthJJbqjSPfPmLmJoQTlGl\nhTBeaZb3hI7MjFioHo991K3TrN5TwJjH1zP6sfXMf/ZLfvTfr/jzhkze21/IgYIKKuuanGCsHUip\nQmJhcTDzIfeM6WS+/PJLRo8ezfjx4zl69Cjp6eldOs+KFSs4ffr80qvly5dz5EjH/U+8AZ1T8BbM\nicqybBg8mbzSGiYkRXvWJhdjrqw6bqwhPUZCdbH3JJnNxA9TNh/7ECbf3aVTfPxNEQ+/e5BJyTGM\n6h/FcWM135yq5ONvirBcxhEXHkhyXJjpJ5yU+DBS4sIYHBtKkL+TVhl/8zYU7ILrnoPgKOec0wVI\nKZFSYjBcfN/75ptv8sADD7Bs2TLy8/O7PMaKFStIT09nwIABALz00ktdPldPQjsFb6FPEhgCwJhF\nXWMLpyrq+HbGoM6P82IshfHSW03hF2+bKYAKIW39q2o+4+BK7C+zS7l31X4uHRzNq0snERJ4/uLe\n2NzKyfJa8kqrOW6s4bixhrzSGj47VoqxurBtP4OAgdEhDO8bwa/mjyK5q2XMDdWw8VGVJxl/a9fO\n4ULy8/OZO3cus2bNYseOHdx33308//zzNDQ0kJqayquvvsqqVatYvXo169evZ9OmTTz55JNtx7e0\ntPDwww+zefNmGhoa+OEPf8g999wDwB//+EfeeOMNDAYDc+fOJSMjg71793LrrbcSEhLCjh07mDt3\nLs888wwZGRmsXLmS3//+90gpmT9/Pk8//TQA4eHh3Hvvvaxdu5aQkBDef/99+vbt65H3yxbaKXgL\nfv4QmwrGnDZJaV+Ut7DkAmE8ae625mUzBVAtXL/8E2R+ApfYfzHdd+Isd7++j7SECF5eOvEChwAQ\n6G8gLSGcNCtlyefqm8g3OYk8k8PYklXKkhd28tY9U0iK7cJn58s/QVURfPt1sHIHbuY3Hx7myOlz\njp+/A0YNiOSxhZ33Z8jMzOTVV1/liSeeYPHixWzatImwsDCefvpp/vznP/Poo4+ydetWFixYwI03\n3njBTOHll18mKiqKPXv20NDQwLRp05g9ezbHjh1jzZo17Nq1i9DQUMrLy4mJieG5555rcwKWnD59\nmoceeoh9+/YRHR3N7NmzWbNmDYsWLaKmpoYpU6bw5JNP8vOf/5wXX3yRX//61059r7qLdgreRGwa\nlGa2SUr7mhBee4ID/EiMNgvjZauy3D5JnjbLcQZcApGJqgrJTqdwtOgcy17dTb+oYF6/axJRIY7V\nnEcGBzA2sQ9jE8/rSx4tOseSF3ey5IWdrLp7KoNjHeidfGof7HgOxi2BQZMcssWdJCUlMWXKFNau\nXcuRI0eYNm0aAI2NjUydOrXDYzds2MDBgwf53//+B0BlZSXZ2dls2rSJZcuWERqq3q+YmI5ne3v2\n7GHmzJmYxTtvvfVWtmzZwqJFiwgMDGTBggUATJgwgY0bN3br73UF2il4E3HDIGsd+WcqAboeBvAi\nUuJNFUgiW7Un9fPCj6wQas3CV6+pEExQx84831jD7S/vJizInze+O4n4iCCnmDGyfyRvLp/MLS/u\nYsmLO1l19xQGxdjhGE7vhze+pUqB7dA3sueO3lWEhanvhJSSa665hpUrV9p9rJSSv//971x77bUX\nvL5u3TqHCjo6qhALCAhoO5efnx/Nzba7KnqKDquPhBBT3GWIxg7ihkJrM5VFWQzsE3JROMEXSTWt\nVZBl2Wqm5K2MuVGp3b7/ww6b7xRX1nPby7tolZI3vjuZxGgH7ubtYPSAKN5cPpmq+iaWvLiTUxV1\nHR9Q9DW8vgiComDpWghPcKo9rmLKlCls27aNnJwcAGpra8nKyurwmGuvvZZ//etfNDWpSq6srCxq\namqYPXs2r7zyCrW1tQCUl5cDEBERQVVV1UXnmTx5Ml988QVGo5GWlhZWrlzJjBlu6u3uBDorSf2n\nW6zQ2Icpni6NWT6fTzCTEh9Gc1MjlB/3znyCmUGTYPaTcGQNrL3PqqT22ZpGbn95FxW1Tby2bJLV\nXIEzSB8YxX+WT6ayroklL+zktC3HUPwNvH49BEXA0g+hz2CX2OMK4uPjWbFiBUuWLGHs2LFMmTKF\nY8eOdXjM8uXLGTVqFJdeeinp6encc889NDc3M2fOHK677joyMjIYP348zzzzDABLly7l+9//PuPH\nj6eu7vx72L9/f5566ilmzZrFuHHjuPTSS7n++vZNJ3suoqOpjhDiKynlpW60p1MyMjLk3r17PW2G\nZ6irgKeTeKb1FqozfsTj1/l+Y/QduWX88qX3+Dzoflj0Lxh/i6dN6h6f/ha+fAam3XtBKKa6oZlb\nX9zJseIqXrtrElNSYl1uyv6TZ7nj5d3Ehgey6u6p9IuyWAh55gi8tgD8g9UMISalw3MdPXqUkSNH\nuthijb1Y+38IIfZJKTNsHNJGZwHaFCFE+x4IbUgpr7O1TeMCQvrQEhrP4HOnaOglM4XUhDBShGmB\nUKwPiMpd+Wul9Lrtb6rB/RU/o76phe+9tpdDp8/xwu0T3OIQAC4ZHM2KuyZxx8u7uMWUY0iIDIaS\no/DaQpXYv/PDTh2CxrfozCmUAn9yhyEa+6gOTyaluohGH688MhMfHsSogDPqSZwX5xTMCAFz/6gc\nw6e/oSUokh8dHc/O42X85dvjuWqke2vWJyQpx3DnK7tZ8uJOVt8QQ+zbi1XTnDvXqjJoTa+iM6dQ\nLaX8wi2WaOziTOAgUsURGnqJUxBCMDa4hMqmPkSF+MgKboMBFv0TWV+J+PgBQhp/yG+u+x6LLvFM\nk5qJQ2J4delEHn/1fXhtKa0h/hiWrvUNJ6xxmM4SzeVusUJjN3lyANGimr7+Vpra+yiphiLy5ABP\nm+FUpMGfJ8MeZnfrCP4a9C/uiO24MsbVTI6q4P3wp6C1hXsMj1MW4oXrQTROoTOn0CiE2CmE+IMQ\nYqYQwguLxH2Lb+pVeEEYsz1sifvo31zA0aa+1Db2vJrurvKXTdm8tLOILRP+jqHfGFh9O5zY7hlj\nyvNgxQICRQsnF65iS0Ust760i/KaRs/Yo/EoHToFKeVcYCawGfgWsFMI8a4Q4m4hhPfUp/kQu6tM\nqynLeolTqC0npKmCPDngQhltL+blrcd59tNsvp2RyIPXZSBue0eVe/73O2pdgDspPw4rFkJzHdzx\nAZdkTOOlOzPIM9Zw20u7qKjVjqG30al0tpSyXkq5Tkp5r6mc6X5ULuI5IcRul1uoaaO2sZl95yJp\nEUoYr1dQphYf5coBbfIe3szbewv47dojzBndj99/a4xa3RoWB7e/p1RH31is2ne6g7MnVJVRYzXc\n8QH0UxLSVwyN58U7Msgpqea2l3dRWesmWe5uUFFRwT//2TOWVT366KNs2rTJoWP27t3LT37yExdZ\n5BgO91OQUh4H/g2sBi53ukUam+SV1tCKgdqIIWDM8bQ57sHk/I7L/l4/U1h3qJiH3jnIFUPj+NuS\n8fj7WXz9ohLhjvdVddLri6CiwLXGVBSodQgN59S4/cdesHnGsHj+ffsEsoqruf2VXe7r19BFbDmF\nlhbbq8ddxRNPPMHVV1/t0DEZGRk8++yzLrLIMTqTuYgUQvxCCPGcEGK2UPwYyANuklLquaUbMd8p\ny7ihvWemYMwGQwCiTxJ5Ru91CtlnqvjJyv2MG9SH52+bYL2/QWwq3PYuNFTBG4ugutQ1xlQWKodQ\nVwm3r4EB463uNmtEAv+89VKOFp3jjld2c66+5zqGhx9+mNzcXMaPH8/EiROZNWsWt9xyC2PGjAFg\n0aJFTJgwgdGjR/PCCy+0HRceHs6vfvUrxo0bx5QpUzhzRpU/v/3226SnpzNu3DimT58OqP4JixYt\nYuHChSQnJ/Pcc8/x5z//mUsuuYQpU6a0yV8sXbq0TVTv4YcfZtSoUYwdO5YHHnjA5rk3b97cJpRX\nXl7OokWL2lZiHzx4EIDHH3+cu+66i5kzZ5KSkuIyJ9JZ4vgN4CywA1gOPAgEAtdLKQ+4xCKNTfJK\naxACQvuPhOProbkR/AM9bZZrKcuB2FSSQiPJLfHe8NGaA6dokZJ/3z6BsKAOvnb9x8Ktq9Vs4T+L\n1WpiZzazOXdahYxqy1XIamDHggVXj+rLP265lB+8+RVLX9nNG9+d3LH9AJ88rCQynEm/MTD3DzY3\n/+EPf+DQoUMcOHCAzZs3M3/+fA4dOtTWvP6VV14hJiaGuro6Jk6cyA033EBsbKxNKesnnniC9evX\nM3DgQCoqKtrGOXToEPv376e+vp60tDSefvpp9u/fz09/+lNef/117rvvvrZ9y8vLee+99zh27BhC\niLbz2Dq3mccee4xLLrmENWvW8Nlnn3HHHXdw4IC63B47dozPP/+cqqoqhg8fzv/93/8REOCYgm5n\ndBY+SpFSLpVS/htYAmQAC7RD8Ay5pdUkRofgnzAMZAucPe5pk1yPMQti00iJC+e4sYbW1h7Qo7gL\nrD98hklDYkiIsKOn9uAp8J3/qJXF/70ZmjoRrbOXqmLlEKpL4LZ3ILFTxQMAZo/ux3O3XMJXJyt4\n9jPvKHCYNGlSm0MAePbZZ9tmAwUFBWRnq7+jvZS1ub/CtGnTWLp0KS+++OIFIahZs2YRERFBfHw8\nUVFRLFyo2q2OGTPmoi5ukZGRBAcHs3z5ct5999026W1b5zazdetWbr/9dgCuvPJKysrKqKxUysjz\n588nKCiIuLg4EhIS2mY2zqSzmULbfFFK2SKEOC6lvFgWUOMW8kprSIkLhzjTnaMx2zs7kdlLS7Oq\njhkxn9SIMOqaWig6V8/APiGetswhckurVdJ2sgMFe0OvhsUvwP/ugtV3ws1vgp+Dd4QNVWqmVZar\nHg+uVo7htncc7okwJ70/N01I5JWtx/lOxqC2rnhW6eCO3l2YJbRBhWY2bdrEjh07CA0NZebMmdTX\n1wO2payff/55du3axUcffcT48ePb7tSDgs7LmBsMhrbnBoPhIhlsf39/du/ezaeffsqqVat47rnn\n+Oyzz2ye24w1PTqzjZbju0p6uzOnME4IYW6hJIAQ03MBSCllpNMt0liltVWSZ6xWujixiepFX88r\nVJyA1iaIG0ZqpKk1Z0m11zmF9YeLAXXH7RDpi1Ui+MN74b3vKydhaJeLaGmCs/mmi3+OulEwO4Hq\nYosdhdIwuvVtSOq42Ywtfj5nBOsOFfPbtUd4dVnParRjS8YaVLOc6OhoQkNDOXbsGDt37uz0fLm5\nuUyePJnJkyfz4YcfUlDgeOK/urqa2tpa5s2bx5QpU0hLS7Pr3NOnT+fNN9/kkUceYfPmzcTFxREZ\n6b5LbYdOQUrp+4L9XkLRuXrqm1qVZHZwJIT3ayvX9FnMTi92KGnRyinklFQzfVi8B41ynPWHihmX\nGMWArjizCUuVOu6mx5RiaeIEddE3Zqv//9l8FUo0Exqr+k6kXaUezT8xyRDQPWcaHxHET64aypMf\nH+WzY2e4ckTP6S0cGxvLtGnTSE9PJyQk5IK+x3PmzOH5559n7NixDB8+nClTOm8T8+CDD5KdnY2U\nkquuuopx48ZddEffGVVVVVx//fXU19cjpeQvf/mLzXN/8cV5NaHHH3+cZcuWMXbsWEJDQ3nttdcc\nGre7dCid3RPprdLZW7JKueOV3az83hSmpsbCigWqactyx+qhvYptz8LGR+ChfGRwHy757Ubmpvfn\nqcVjPG2Z3RRV1jH1qc948Nrh/HBWN7SENj0OW9VFBf8Q08U+9fxFP26omgmEdtwqsrs0Nrcy529b\nkBLW3XdFWxWVls7uWbhSOlvTQzhapKJ4bY1X4obCoXdUsxYHWgV6FWXZEBYPIdEIIC0+3OsqkDYc\nVonAOekOho7ac9VjMPY7quFNxAAlqucBAv0NPLpgFEtf3cOr2/L5/gytoupreOaTpXGYHXllpMSH\nne/XGztUyS/XGD1rmCsxZl/QQyEtIZwcL1vVvO5QMWkJ4aR2V9VWCEgYqRa5ecghmJk5PIGrR/bl\n759mU3Ku3qO2aJyPdgpeQFNLK7uPlzMtNe78i+bWlL6cbDZmXyDfnJYQTnlNo9cItZ2taWR3fjnX\nju45sXdn8ciCkTS1SP6wruMWlxrvQzsFL+DrggpqG1u4LNWiI5f5Yumrwni15VBrvKAvszl0luMl\nIaRNR8/Q0iqZM7q/p01xOkmxYSy/Ipl3vzrFvhNnAeullBr3093/g3YKXsD23DKAC9s0Rg1S1Si+\nKqFtrqxqFz4CyC7xjqUy6w8XM7BPCOkDfbNy+4ez0ugbGcTjHxwmKCiYsrIy7Rg8jJSSsrIygoPt\nWCRpA51o9gK25xoZ1T+S6DALSQuDH8Sk+q5TMP9dceedwoCoEEIC/LxiplDT0MyWbCO3TBrctvDI\n1wgL8ueX80Zy76oDfFWWyARRRWmpi/SaNHYTHBxMYmJil4/XTqGHU9/UwlcnKrjzMiudsOKGul9/\n312UKSE8+pz/uw0GoZLNXuAUvsgqpbG5lWsdXbDmZVw3bgBv7DjBU+uy+eyBmUSFOFeHR+N+XBo+\nEkLMEUJkCiFyhBAPd7DfjUIIKYSwT4ylF7HvxFkaW1q5zDLJbCZumFr129zgfsNcjTFb1d37XXjf\nkpbgHWWp6w8XExMWyMQhPtJX2gZCCB6/bjTltY08+6mPzlp7GS5zCkIIP+AfwFxgFLBECDHKyn4R\nwE+AXa6yxZvZlmPEzyCYmGxlUVLcUJCtqp2ir2HMviB0ZCYtIZzTlfVUN/Tc1pyNza18drSEq0cm\nXNgzwUdJHxjFzRMH89r2fLLPeEe+R2MbV35iJwE5Uso8U9+FVcD1Vvb7LfBHQBc8W2F7bhnjEqMI\ntyZXbL5o+lpeoaVZOTorTsFc79+TZwvbc41UNTT7fOjIkgdmDyM00I/ffHhEJ5u9HFc6hYGApdJT\noem1NoQQlwCDpJRrOzqRqSf0XiHE3t6UyDpX38TBwgrroSNQ8gbge2sVzEJ4sdZnCtCzy1LXHz5D\nWKAf09Js/N98kNjwIH56zTC25hjZcMT5cs4a9+FKp2Ct5KLtFkIIYQD+gur53CFSyheklBlSyoz4\neO8SQ+sOe46X0yrhsrRY6zuYJQ98TRjPSuWRmaTYUAL8RI9d2dzSKtl4pJiZIxIIDuhdepK3TUli\nWN9wfvfREeqb3N8GU+McXOkUCoFBFs8TgdMWzyOAdGCzECIfmAJ8oJPN59meW0agv4FLB3eQrIxL\n872ZgnlBXuzFAnIBfhG1TdgAACAASURBVAaGxIb12JnCVyfPYqxu7FWhIzMBfgYeWziagvI6XvrS\nB/NcvQRXOoU9wFAhRLIQIhC4GfjAvFFKWSmljJNSDpFSDgF2AtdJKXukBOpxYw1/3ZTl1s5f23KM\nZCRFd3zHGTcMjDlKGM9XMGZBaJxNxc+eXJa6/lAxgX4GZg3vPTNaS6alxTE3vR//+DyX0xVO6hin\ncSsucwpSymbgR8B64CiwWkp5WAjxhBDiOleN6yr++XkOf92UzVcnz7plvLLqBo4VV10obWGN2KHQ\nUKlaLPoKxhyroSMzaQnhnCiroaG5Z4UopJSsP1LMtLRYIoJ7b73+L+eNpFVKnvpE6yJ5Iy6tl5NS\nfiylHCalTJVSPml67VEp5QdW9p3ZU2cJDc0tbd2zPvqmyC1j7swrB2CqrSSzGfPF05c0kMqsl6Oa\nSUsIp1VCvrHWjUZ1zpGicxSU1/XK0JElg2JCuWdGKh9+fZrdx8s9bY7GQXy/iNoJfJll5Fx9M/ER\nQXzyTbFbQkjbc42EB/kzLjGq4x3bylJ9JK9QdxZqSq1WHpnpqRVI6w+fwSDg6lG+p4rqKP83I5UB\nUcE89sFhWtwYctV0H+0U7GDtwdP0CQ3gwWuHU3yunv0Frg8h7cgtY1JyTOeLnyITVScuo49UIJn/\njg5mCqnx4QjR85zChsPFZAyJIS48qPOdfZyQQD9+NX8UR4vOsXL3SU+bo3EA7RQ6ob6phY1HzjBn\ndD/mpPcj0M/ARweLOz+wGxRV1pFnrOk8nwCq4UqsD1UgmcNgFpLZ7QkO8CMxOqRHqaXmG2s4VlzV\n60NHlswb048pKTE8syGTilrv6IGh0U6hUz4/VkJNYwsLxw0gMjiA6cPi+eRQkUtDSNtzlFT2VHuc\nAqi7al9xCsasi4TwrJEW37MqkMw5J19sqNNVhBA8tnA05+qa+MtGH/l89gK0U+iEtQeLiAsPZLJJ\ne2j+2H4UVdazv6DCZWNuzy0jOjSAkf3s1OGPGwYVJ6HJB5RCjNkQk3yREF57hvaNIM9Y02Pi1esP\nF5M+MJLE6FBPm9KjGNk/ktumJPHGzhMcKz7naXM0dqCdQgfUNDTz6bEzzBvTvy22f9XIvgT6GfjY\nRVVIUkp25BqZmhqLwWCnDn/cUEBCea5LbHIrxuwOQ0dm0uLDaWxupfCs5yuQzpyr56uTFVw7SoeO\nrPGza4YRGRLA4x8c1rpIXoB2Ch2w6egZ6ptaWTB2QNtrKoQUxyffuCaEdKKsltOV9Z2XolriK8J4\nZiE8KyuZ25PagyqQzFo/c9K1U7BGn9BAHpg9nJ155WzUukg9Hu0UOuDDr4voFxlMRtKFMhPzxvTn\ntItCSNtyjQD2JZnNtAnjeblTMAvhdVB5ZOZ8a84e4BQOF5MSF9Zmk+Zibp44iKTYUP7+WY6eLfRw\ntFOwQWVdE1uySpk/tv9FYZyrR7kuhLQ9t4y+kUGkxIXZf1BgmCpN9fYFbMbOK4/MRIUEEB8R5PGZ\nQmVtEztyy5g9up/Ptt10Bv5+Bn4wM5VvTlXyRVbvUTr2RrRTsMGGw8U0trSyYGz/i7ZFBgdwxVDn\nh5BaWyU7c8uYlhrn+AXGF4TxOhDCs8bQHqCB9OmxMzS3Sl11ZAffuiSRgX1C9Gyhh6Odgg3WHiwi\nMTqE8YP6WN1uDiEdKHReCCmrpIqymkb7S1Et8QVhPGM2hMbaFMJrj7k1pycvMOsPF9MvMphxidY/\nJ5rzBPob+P6MFPadOMuOvDJPm6OxgXYKViivaWRrjpEFYwfYvGO/elRfAvwEHx90Xghpm6PrEyyJ\nHQqNVVDl2oV1LsXOyiMzaQnhVDU0U1LlmR7VdY0tfJFVyuzRfe2vFOvl3JQxiPiIIJ77zEdW4Psg\n2ilYYd2hYlpaJQvHXRw6MhMVEsAVQ+P55FCx0+5Ud+QaSYoN7Vqtuy8I45Vl2x06AlWWCpB9xjMh\npC+ySqlvatWrmB0gOMCPe6ansD23jH0ntFheT0Q7BSusPXialLgwRvXvePHYvDH9OVVRxwEnVCE1\nt7SyK6/cduvNzvB2YTyzEJ6DMwWAHA/JXWw4XEyf0AAmJdsX7tIobpk8mOjQAD1b6KFop9COkqp6\nduaVsWCc7dCRmWvMISQnVCEdOn2OqoZmx0pRLYkYAAFh3iuMZ4cQXnviI4KIDPb3SGvOppZWNh09\nw1Uj+hLQmWih5gJCA/1ZfkUKn2eWcuhUpafN0bRDf5rb8ck3xbRKWGil6qg9USEBXJ4Wx8ffdD+E\ntN20PmFKShedgsEAsaneO1Noqzyy3ykIITzWhW1nXhnn6pt11VEXuWNqEpHB/nq20APRTqEdH359\nmuF9IxjaN8Ku/c0hpK8Lu3fHsz2njOF9I4iP6Ibsctww780pGLOVEF50x0J47fGUU1h/uJiQAD+m\nD+udbTe7S0RwAEunJbPucDGZxT1H7VajncIFnK6oY++Jsx0mmNsze1S/boeQGppb2JNf3rWqI0vi\nhkFFATR5YW9cY5ZJCM+xNpZpCeEYqxvdKs3c2irZcPgMM4fHd9w/W9Mhyy4bQligH//4XM8WehLa\nKVjwkam81FLrqDOiQgOYlhbHRweLuhxC2n+ygobmVqaldTHJbCYuDZBQ5oXCeGU5DoWOzHiiC9uB\nwgpKqhp01VE3iQ4L5LapSaw9eJo8D+SFNNbRTsGCtQdPM2ZgFEMckZjgfAjpYBdDSNtzyzAIul/F\nYq7c8ba8glkIz4Eks5mhCSrM506nsP5QMQF+glkjEtw2pq+y/PIUAvwM/GuzF97I+CjaKZg4UVbD\n14WVVmUtOuPaboaQduQaGTMwiqgQx0InFxGTqh7LvGw6XnECWhq75BQG9gkhOMDgNqcgpWT94WKm\npsZ1//+lIT4iiCWTBvPe/lMUlHteBl2jnUIba02ho/ldcAptIaRvHA8h1TQ0s/9khWNS2bYIDIWo\nwd43UzA7sS6EjwwGQUpcuNvUUrPOVJNfVqurjpzIPTNSMAjBv7fo2UJPQDsFEx9+fZpLB/fpcues\neWP6U3i2jm8crLvek19Oc6tkWlo3k8xmvFEYr00d1XGnAO6tQFp3qBgh1BoVjXPoHxXCjRmJrN5T\nSHGlD3QP9HK0U0DFo48VV/3/9s48PKoqzf+fNwlJSEL2hRD2JAJhDwgouKCIgiBqa6tjK/q4tGv3\nPDM90073PNptPz1jt7/5td3t3uq4NGp32+2GCwYFkV3WACGQhB1CNkLIQtY688e9FYsiS1Xq3qok\ndT7Pkyc35557z1u3buqt+573fF+vJpjdmZ+TRliI8ImXIaQNJVUMCBWmj7BoVWxfFMar3O+VEJ47\n2akxHD99lobmVosNO58Ve04ybXgCqYMibR8rmHjwskzalOLlNQcCbUrQo50CxgSzSM9CR07io8KZ\nnZXMp16GkNaXVDF1eAIDwy1KbUzKgpZ6OHPCmvP5gx5mHjlxZiAdqKi3yqIOOXqqgYLSMzrryAaG\nJUZx/ZQM3t58mMq6wAgcagyC3ikopfh45wlmjEwkLda3b3/XTkzn6Kmz7D7uWYHymoYWdp+o6bm0\nRUc4M5DsWMTWatM/a2VRj0NH4L+01BV7DAVa7RTs4aG5mTS1Onh17cFAmxLUBL1TKDxZS0lFPYsn\n9zx05GT+eO9CSBsOVKEUPRfB6wi76jXveR+eGgGnLP6HPXsa6st9cgojkqIJDRGKbBbGW7HnJOPS\nYxme1LN5J03XZKbEsGjSEN5cf8ivixE15xL0TuHjnScIDREWWFB0PT4qnIu9CCFtKKlk4IDQTgv5\n9IhB6RAeY61TUAq++f/Qeha2vWHdecGnzCMn4WEhjEiKsvVJoaK2iS2Hq3XWkc08PDeT+uY2Xl9/\nKNCmBC1B7RSUUizPL+XizCSSYnzQHHLh2omDOXKqgT0nug8hrS+p4sJRiYSHWfg2iBjzClZmIB3Z\nACfzISIWdrwNbS3WnduLusxdYXdpzi/3lqGUDh3ZzdjBsczPSeO1tQepbbTwPtN4TFA7hfxjNRw5\n1cBiH7KO3JmfM5hQD0JI5bWNFJXXWTuf4CT5AmsXsG18ASLjYfEzUFcG+1dYd+7K/RAS5rUQnjtZ\nqTEcrmqgudVhkWHnkldQxtCEgYwd7JlQoqbnPHJFFmcaW3lr4+FAmxKUBLVTWJ5/ggGhYum3v4To\ncC7OTOo2hLShxCi9aZtTqDkKzRZk41QfhsLlMG0pjFtihKe2ven7eZ1UFUGC90J47mSlxtDqUByu\nsj4DqaG5lbXFlVyVk9ZtjQ2N70waGs9lF6TwyjcH/ZJmrDmXoHUKDocROro0O4W4KGvlCq6dmM7h\nqq5DSOuLq4iNDGP8kDhLxwZMYTysEcb79k+AwIX3QWgYTLkdivOg5rjv5wZjTYWPoSOArBT7NJC+\nKaqkqdXBVeP0fIK/+NGVWZyqb+adzUcDbUrQEbROYduRakprGlnkhUy2p8wf330Iaf2BSmaNTiLU\njoLvVgnjNdUZTwXjFkP8MKMt9w5QDtixzLdzAzja4FTJd07MBzJTDRFDO5zCyoIyYiPDuFCX3fQb\n00YkctHoJF5eU0JjS1ugzQkqgtYpLM8vJSIshHk2fPtL7CaEdPRUA0dPnbUndASQOBoQ3+cVdr4D\njTUw68Hv2hJGwui5sO0tcPgYv3cK4fmQeeQkKjyMjPiBlpfmbHMoviosZ+7YVF120888ekUWZWea\neG/rsUCbElTYepeLyDUisk9EikXksQ72/4uIFIhIvoh8KSK+zTZ6SJtD8cmuUuaOSWVQpD1Klwu7\nCCG1zyf4Wj+hMwYMhHgfhfEcDtj0EgyZCsNmnrsv906oOQIHVvlmZ3tdZt/DR2CPBtL2I9VU1Tdr\nraMAcFFmErnD43lhdQktbfYkEGjOxzanICKhwHPAAiAHuE1Ecty6bQemK6UmAe8Bv7XLHlc2Hayi\norbJkgVrnXG1GULqSE57XUklyTERZJsrcW0hOdu3tQolXxmTwDMfNNJcXRl7raFV5OuaBafT8mHh\nmitZqTGUVNThcFin+5RXUMaAUOEyXXbT74gIj16RzfHTZ3l/u0VzWJpusfNJYQZQrJQ6oJRqBt4F\nlrh2UEqtUko5RdQ3AkNttKedj3eWEhUeyhU2FklJjA7notHnh5CUUqwvqeLizCR7M1mcaak9DfFs\nfB5i0mD8DefvC4uAybdB4adQV9FzG6uKYGBij4Xw3MlKjaGxxcHx09aVI80rKGPW6CTbnig1XXP5\nmBQmZMTy/Kpi2ix09prOsdMpZACuqQPHzLbOuAf4rKMdInK/iGwRkS0VFT58CAEtbQ4+313KlePS\nrBOh64RrJ6VzqMoQUXNSUlFHRW2TffMJTpKyoKUBansgjFexH0q+hAvvhbDwjvvk3gmOFmPeoadY\nlHnkxGoNpJKKOg5U1jNfh44ChojwyNxsDlU1sDy/D4k89mHsdAodfQ3u0NWLyA+A6cDTHe1XSr2s\nlJqulJqekuLbY/y64kqqG1pY7IMiqqd0FEJa374+wab5BCe+ZCBtehFCw2Ha3Z33SRkDw2YZ2Uk9\nlemu3G9J5pGTrBRrnUJeQRkAV+pU1IAyPyeNC9JieG5VsaWhQU3H2OkUjgHDXP4eCpzn6kVkHvBz\n4DqllO2aucvzSxkUEcZlY+yPEX8XQjrZHkJaX1xFRvxAhiUOtHfwdmE8LzOQzlYb3/4n3gwx3Vyj\naUuNENCRDd7b5xTCsyDzyElCdDjJMeGWOYWVBWVMyIhlSLzN75WmS0JChIfnZrG/rM7reiUa77HT\nKXwLZIvIKBEJB24FPnLtICJTgZcwHEK5jbYA0NTaxoo9J5k/fjARYfaGjpwsnJjOwcp69pbW0uZQ\nbDhQxewsm+cTwJgPiIj1/klh25tG2GnmA933zVlijLG1BxPOVdZmHjnJTImxRC21sq6JrUeqbUlZ\n1njPoklDGJcey1OfFep1CzZjm1NQSrUCjwArgL3AX5VSe0TkSRG5zuz2NBAD/E1EdojIR52czhLW\n7K+ktrHVlgVrnXH1+LT2ENLe0jPUnG2xP3QEPRPGa2uFzX+CEXMgfVL3/cOjjSeKgg+Mb/7e4GMJ\nzs5wpqV6Wyvbna8Ky1FKl93sLYSGCI8vyuH46bP8SVdnsxVb1ykopT5VSl2glMpUSv3abHtcKfWR\nuT1PKZWmlJpi/lzX9Rl9Y3n+CeKjBjDHrvUBHZAUE8Gs0Yl8uquUdcWVgJF/7Re8FcYrXG5oJs3y\n4CnBSe6d0NoIu/7mnW1VRaYQ3kjvjuuG7NQYzjS2UuFj9a68gjIy4geSkx5rkWUaX7koM4lrxg/m\n+dUlupazjQTNEs2zzW3kFZSxYMJgv69MXTgxnQOV9SzbdITMlGifK7x5THI2nDluyFV4wqYXjUVv\nYxZ6PsaQKZA+2QghefPtvHK/JUJ47mSl+q6BdLa5jW+KKpg3LlUL4PUyfrZwHG1K8dvPCwNtSr8l\naJzCqn3lNDS3schCmWxPuXr8YEIEjpxq8E/oyIkzNOPJ08KJ7caE8YwfQoiX8y25d0LZLuMcnlJZ\nbHnoCKxJS11XXElji4N5OnTU6xieFMW9c0bxj+3H2X6kOtDm9EuCxim0tDmYMiyemQEQNUuOiWDW\naCNkNDvLT6EjcElL9WBl88YXYUA0TP2B9+NMvBnCBnouqe0UwkuyLh3VSVpsBDERYT45hbyCMgZF\nhDFzlB/fK43HPDQ3i5RBETy5vMDnuSPN+QSNU1gyJYMPHp5NWIBEzW65cBixkWHtzsEvJI4GCTHi\n911RWwa7/w5Tb4eBPSgNGhlnrHze9Z5noSqnEJ7FmUdgLHbK9EEDyeFQfFlYxuVjU62tiKexjJiI\nMP796jFsP3KaD3foBW1Wo+96P7FkSgbbH59PfFQnK4TtICwC4kd0n4G05TVjdfKMH/Z8rGlLobkW\n9rzffd92ITzrw0fgW2nO7UdPU1nXzLxx9kmgaHzne7lDmTQ0jqc+K9SFeCxGOwU/YkvthO5Izu56\nAVtrE2x5FbLn+7a6eNhMSB7jWQipypq6zJ2RlRpDeW0TNWe9r/G7cm8ZYSHC5WO0U+jNhJgpqifP\nNPLi1zpF1Uq0U+jvdCeMt/vvUF9xbs2EniBiTDgf2wzle7vuW7nfUiE8d3yRu8grKGPm6ETiBmoB\nvN7O9JGJLJ48hJe+LuFYdUP3B2g8QjuF/k5SFrSehTMdFCpRCja+ACljjcI5vjL5NggZ0P3Tgk2Z\nR06cGUglXjqFg5X1FJfX6bKbfYjHFoxFBJ76TKeoWoV2Cv2droTxDq+Hk/kw84fn10zoCdFJMG6R\noZ3U0sXioqoiW53CsMQowsNCvK7CttIUwNOpqH2HjPiB3H9pJsvzS/n20KlAm9Mv0E6hv9OVMN6m\nFyAyHibdat14uUsNUb3C5R3vb6yBujJLhfDcCQ0RRidHex0+yisoY1x6LEMTomyyTGMHD1w2msGx\nkTz5cYFWUbUA7RT6O9EpRsqo+5NC9WEo/ASm3QXhFn4IjrrMyHjqrCqbzZlHTrwtzXmqvpkth09x\nlc466nNEhYfx2IKx7Dpew3vbdD1nX9FOob8jYnwrd1+rsPllQGDGfdaOFxICuXfAwTVwqoOsEJsz\nj5xkpcZwtLrBY0XNrwrLcSi4KmewrXZp7GHJlCFMHR7P0yv2UdekU1R9QTuFYCD5gnNXNTfVwba3\nIOc6iLOhAuqU241Fc9veOn9f5X5bhPDcyUqNQSmjeponrCwoY3BsJBMytABeX0REeGLxeCpqm3hu\nlZc1RDTnoJ1CMJCcBbWl0GTWGdj5DjTVwEwf01A7I3YIZF8NO5ZBm9tagcoiwyFYLITnjjcaSI0t\nbawpqmBejhbA68tMGRbPjbkZvPrNQY5U6RTVnqKdQjDgqoHkcBhqqEOmwrAZ9o2Ze6cxoVz0xbnt\nVdbWZe6MUcnRhIhnaanrSyppaG7ToaN+wE+vGUtYqPBfn3azVkbTKdopBAPOD+GqYij50vg96yFr\n0lA7I3s+xAw+tyqbow2q7BHCcyciLJQRSdEUeeAU8grKiYkIY9Zo/4slaqwlLTaShy7P5PM9J1lf\nUhloc/ok2ikEAwmjQEKNeP7GF4wP65zr7R0zNMwQ2CvOg5rjRtvpI9DWZHvmkZPMlO4zkBwOxcq9\nZVx2QYrfSrRq7OXeS0aTET+QJz8uoE2nqHqNdgrBQFi4Ecff97nxpHDhPUab3Uy9A5TDmFsA2+oy\nd0ZWagyHquppbetE4gPIP15DRW0T83J0Kmp/IXJAKD9bOI7Ck7W8++2RQJvT59BOIVhIzjYK4YRG\nwLS7/TNm4ihj3cK2t4y5DOdaCRsXrrmSnRpDS5vi8KnOJx3zCk4SGiLM1QJ4/YqFEwczY2Qi//PF\n/h4JIwYz2ikEC86QzcSbISbFf+NOWwo1R+DAKmOie2CCIYfhBzzJQFpZUM6FIxP8K2musR0R4fHF\nOVQ3NPPHLz0oMqVpRzuFYCF9irF2YNYD/h137CJDEXXbm37LPHKS2Y1TOFLVwL6yWp111E+ZkBHH\n96cN4/X1hzjgpQ5Wb+ToqQa/VJrTTiFYGH8D/HgnDJ7o33HDIgz11MJPDPE9P4WOwKjQlR4X2alT\nyNtrCOBpVdT+y0+uHkPkgFB+/UnfTlH9cMdx5v9uDW9uOGz7WNopBAshoRA/PDBj595pVHZrrPGt\nkE8P6EoDKa/gJGPSBjE8SQvg9VdSBkXwyBVZfFlYzpr9FYE2x2uaWx088eFufvzuDiZmxLFggv1P\ntdopaOwndaxRmQ38Gj4CwymUVNSdp555uqGZbw9Vc5WWye733D17JCOSovjV8gJaushE622U1pzl\nlpc38MaGw9w7ZxTL7ptJamyk7eNqp6DxDzPuN9ZKpE3w67BZqTE0NLdReubc+g6r9pXT5lC6dkIQ\nEBEWys8XjqOovI4fvbO9TwjmrSuu5No/rGX/yVqe+6dc/nNRDgNC/fNxrZ2Cxj9MvAn+tRASRvh1\nWGdpzqKy2nPa8wrKSB0UwaSMOL/aowkMV+Wk8fOF41ix5yQ3PLeOg5X1gTapQxwOxXOrirnj1U0k\nRofz4SNzuHZSul9t0E5B4z9i/L8WoKO01KbWNr7eV8GV49IICdECeMGAiHDfpaN5656ZVNY1cd2z\na/nSTDToLdScbeH+t7by9Ip9LJyYzocPz26/f/2Jdgqafk1STAQJUQPOkdDeUFJFfXMb83XoKOiY\nnZXMx4/OYXhiFPe8sYVnVu7vFdXa9pyo4bpn17J6XzlPLM7hj7dNJToiLCC2aKeg6fdkpw4650lh\n5d4yosJDuSjTP4voNL2LoQlR/P3Bi7kxN4NnVhZx/1tbONMYuFXP7209xo3Pr6expY1375/F3bNH\nBVTCXTsFTb8nMzWGovI6lFIopVhZUM6l2SlEDtACeMFK5IBQ/ufmyfxicQ6r91Vw/bPrzpt3spvG\nljb+4x+7+MnfdjJ1eDzLH72E6SMDr9SrnYKm35OVGsPphhaq6pvZdbyGk2cadSqqBhHhrtmjePu+\nWZxpbOH659bx2a5Sv4x9rLqB77+0gXc2H+GByzL58z0zSRkU4Zexu0M7BU2/x3WyeWVBGSECc8dq\nATyNwYxRiSx/9BKy0wbx4LJt/PbzQlslt7/eX8GiP67lYEU9L90xjccWjCXMT+mmntB7LNFobCLb\nxSl8UVDG9JGJJEZrATzNdwyOi+QvP5zFbTOG8fzqEu5+/VtONzRbOobDofj9yiLu+t/NpA2K5KNH\n53D1+N6nuxWY6W2Nxo+kx0USHR7K6n0VFJ6s5ecLxwXaJE0vJCIslP++cRKThsbzxId7WPzsWl76\nwXRyhsT2+JyNLW3sLT3D7uM1fLb7JOtLqrhhaga/vmECUeG98+O3d1ql0ViIiJCZGsNKMy9dr2LW\ndMVtM4YzZvAgHvzzVm58YR2/+d4klkzJ6Pa4huZW9paeYdexGnafMBxBUXldeygqMTqcXy0Zzw9m\njQhodlF32OoUROQa4PdAKPCKUuopt/0RwJvANKAKuEUpdchOmzTBSVZKDPnHashKjWFUcnSgzdH0\ncnKHJ/Dxo3N4ZNl2fvzuDnYdqzkn9l/X1ErBiTPsOl7DbvOnpKIO51REckw4EzLimDcujQkZcUwc\nGseQuMhe7Qyc2OYURCQUeA64CjgGfCsiHymlCly63QNUK6WyRORW4DfALXbZpAlenLUVdNaRxlNS\nB0Wy7L6Z/PqTvbyy9iA7j50mPW4gu0/UcLCyHmdpg9RBEYaC6cR0JmbEMTEjjrTYiD7hADrCzieF\nGUCxUuoAgIi8CywBXJ3CEuAX5vZ7wLMiIsoflSQ0QcXkofGI4BfpYU3/YUBoCL+4bjwTM+J44qM9\nHKs+y4SMOJZMzmDi0FgmDInzi3KpP7HTKWQAR13+PgbM7KyPUqpVRGqAJKDStZOI3A/cb/5ZJyL7\nemhTsvu5ezF9xdY+Zefkp7rtF2j6yvWEvmOrpXZutOpE52P39fRIjdJOp9DRs5P7E4AnfVBKvQy8\n7LNBIluUUtN9PY8/6Cu2ajutpa/YCX3HVm2nd9i5TuEYMMzl76HAic76iEgYEAecstEmjUaj0XSB\nnU7hWyBbREaJSDhwK/CRW5+PgKXm9k3AV3o+QaPRaAKHbeEjc47gEWAFRkrqa0qpPSLyJLBFKfUR\n8CrwlogUYzwh3GqXPSY+h6D8SF+xVdtpLX3FTug7tmo7vUD0F3ONRqPRONHaRxqNRqNpRzsFjUaj\n0bTTb5yCiFwjIvtEpFhEHutg/10iUiEiO8yfe132LRWRIvNnqfuxfrbzdy427heR0y772lz2uU/a\nW23nayJSLiK7O9kvIvIH83Xki0iuyz5/Xs/u7LzdtC9fRNaLyGSXfYdEZJd5PbcE2M7LRaTG5f19\n3GVfl/dMAGz9Nxc7d5v3ZaK5z5/XdJiIrBKRvSKyR0R+3EGfgN+nHtrZK+5TgPZqVH35B2MiuwQY\nDYQDO4Ectz53VmuE9QAACDhJREFUAc92cGwicMD8nWBuJwTKTrf+j2JM0Dv/rvPjNb0UyAV2d7J/\nIfAZxlqTWcAmf19PD+282Dk+sMBpp/n3ISC5l1zPy4Hlvt4z/rDVre9ijKzBQFzTdCDX3B4E7O/g\n/z7g96mHdvaK+1Qp1W+eFNolNZRSzYBTUsMTrgbylFKnlFLVQB5wTS+x8zbgHZts6RKl1Bq6XjOy\nBHhTGWwE4kUkHf9ez27tVEqtN+0AYzHqULts6QoPrmdn+HJv9wgvbQ3kPVqqlNpmbtcCezFUElwJ\n+H3qiZ295T6F/hM+6khSoyOt2++Zj2fviYhzYZ2nx1qBx2OJyAhgFPCVS3OkiGwRkY0icr1NNnpK\nZ6/Fn9fTW+7B+NboRAFfiMhWMaRUAs1FIrJTRD4TkfFmW6+9niIShfFB+neX5oBcUxEZCUwFNrnt\n6lX3aRd2uhLQ+7S/1FPwRC7jY+AdpVSTiDwAvAFc4eGxVuHNWLcC7yml2lzahiulTojIaOArEdml\nlCqx3ErP6Oy1+PN6eoyIzMX4Z5vj0jzbvJ6pQJ6IFJrfkgPBNmCEUqpORBYCHwDZ9NLrabIYWKeU\ncn2q8Ps1FZEYDMf0z0qpM+67OzgkIPdpN3Y6+wT8Pu0vTwrdSmoopaqUUk3mn3/CqOHg0bH+tNOF\nW3F7LFdKnTB/HwBWY3zjCBSdvRZ/Xk+PEJFJwCvAEqVUlbPd5XqWA+9jhGoCglLqjFKqztz+FBgg\nIsn0wuvpQlf3qF+uqYgMwPigXaaU+kcHXXrFfeqBnb3nPvXX5IWdPxhPPAcwwi3Oybjxbn3SXbZv\nADaq7yacDmJMNiWY24mBstPsNwZjcklc2hKACHM7GSjC/gnHkXQ+MXot507gbfb39fTQzuFAMXCx\nW3s0MMhlez1wTQDtHOx8vzH+6Y+Y19aje8aftpr7nTpl0YG6pub1eRN4pos+Ab9PPbSz19yn/SJ8\npDyT1PiRiFwHtGLczHeZx54SkV9haDUBPKnOfRz2t51gTN69q8w7wWQc8JKIODCe8J5S5xYsshQR\neQcjIyZZRI4BTwADzNfxIvApRmZHMdAA3G3u89v19NDOxzHk2J8Xo+hJqzKUKNOA9822MOBtpdTn\nAbTzJuBBEWkFzgK3mu9/h/eMXXZ6aCsYX6y+UErVuxzq12sKzAbuAHaJyA6z7WcYH7C96T71xM5e\ncZ+ClrnQaDQajQv9ZU5Bo9FoNBagnYJGo9Fo2tFOQaPRaDTtaKeg0Wg0mna0U9BoNBpNO9opaPo1\nIrJaRCwthi4i8SLykMvfl4vIcg+PfUZELvVirCEi8l5P7PQFEQkXkTVi1E7XBBHaKWg03hMPPNRt\nLzdMeelZyguJAqXUCaXUTV6MEeqtXZ2M2wx8Cdxixfk0fQftFDR+R0T+XUR+ZG7/TkS+MrevFJE/\nm9svmOJ/e0Tkl2bbAhH5q8t5LheRj83t+SKyQUS2icjfTJ0Z93E77GPq1f/SbN8lImPN9hQRyTPb\nXxKRw6b0xFNApqlv/7R5+hhTaLFQRJaJudrIjZuA9oVH5rj/Zdq0RURyRWSFiJSY+lyIyEgx6xqI\nSKiI/D/TxnwRedTlPI+LyFrgZhGZIoZoYr6IvC8iCWa/1SLyGxHZLEatjkvM9vFm2w7zmGzTxA+A\n2719fzV9G+0UNIFgDXCJuT0d4wN1AIYI2Ddm+8/NFZ2TgMtMXZg8YJaIRJt9bgH+Yn5Q/ycwTymV\nC2wB/sV1QA/6VJrtLwA/MduewKgVkIuhOTPcbH8MKFFKTVFK/ZvZNhX4ZyAHo/bB7A5e92xgq1vb\nUaXURebrfh3DccwCnuzg+Psx5C6mKqUmActc9jUqpeYopd7FkFT4qdlnl/k6nIQppWaYtjrbHwB+\nr5SagvF+HDPbdwMXdmCHph+jnYImEGwFponIIKAJ2IDxYXQJ3zmF74vINmA7MB5D56kV45v2YjPW\nfS3wIcaHaA6wzpQRWAqMcBuzuz5OkbKtGLo/YDipdwFMaYFqOmezUuqYUsoB7HA5hyvpQIVbm1Pa\nZBdGYZVapVQF0Cgi8W595wEvmtcBN1mGvwCISBwQr5T62mx/A6NoTlevcwPwMxH5KYZS61nz/G1A\ns/k+aYIEPYmk8TtKqRYROYShQ7MeyAfmApnAXhEZhfFt/UKlVLWIvA5Emof/BXgYQ7/qW6VUrRmq\nyVNK3dbFsN31cSrotvHd/0VHIaDOaHLZdj2HK2f57nW4H+dwO4ejg3MIncs713fS7s55r1Mp9baI\nbMJwsitE5F6llLOORwTQ6OG5Nf0A/aSgCRRrMD7412A8HTwA7DBF4GIxPuRqRCQNozyhk9UYpSLv\nw/x2jFGparaIZIFR/EVELnAbz5M+7qwFvm/2n4+hpglQi1FW0Vv2Alk9OM7JF8ADzowgc+L6HJRS\nNUC1c74AQ4jta/d+rohRn+OAUuoPGE8uk8z2JKBCKdXig82aPoZ2CppA8Q1GOGWDUqoM49voNwBK\nqZ0YYaM9wGvAOudBZkhjOYajWG62VWCo3r4jIvkYDmCs62Ce9OmAXwLzzTDWAqAUqFWG1v06MYrW\nP93lGc7lEwz10Z7yCoakdr6I7AT+qZN+S4Gnzdc5hY7nJ1y5BdhthtXGYsxJgPH09qkP9mr6IFol\nVaPpBBGJANpMyfOLgBfMyVhfzrkWWKSUOm2JkTYiIv8A/kMptS/Qtmj8h55T0Gg6ZzjwVxEJAZox\nQla+8q/meXu1UxCRcOAD7RCCD/2koNFoNJp29JyCRqPRaNrRTkGj0Wg07WinoNFoNJp2tFPQaDQa\nTTvaKWg0Go2mnf8DgmctqByCzfMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure();\n",
    "plt.plot(wavelength_scan, np.abs(spectra), '-', markersize = 10);\n",
    "plt.plot(wavelength_scan, np.abs(spectra_T), '-', markersize = 10)\n",
    "plt.xlabel('wavelength (microns)')\n",
    "plt.ylabel('R/T')\n",
    "plt.title('sample RCWA spectra for a 1D grating')\n",
    "plt.legend(['reflection', 'transmission'])\n",
    "# plt.axhline(((3.48-1)/(3.48+1))**2,xmin=0, xmax = max(wavelength_scan))\n",
    "# plt.axhline(((3.48-1)/(3.48+1)),xmin=0, xmax = max(wavelength_scan), color='r')\n",
    "plt.ylim([0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1075.3145+8.3458e-14j 1075.305 +5.4572e-14j 1075.3145+2.6694e-14j\n",
      " 1075.305 +1.7704e-14j  869.9639-2.9032e-14j  870.3438-1.1017e-13j\n",
      "  869.9639+3.6598e-13j  870.3438+4.6165e-13j  686.7352-8.5336e-14j\n",
      "  686.5023-1.5441e-13j  524.8468-1.1809e-14j  524.5176-1.4140e-13j\n",
      "  686.7352-2.9469e-14j  686.5023-1.6492e-13j  384.0719-7.4836e-14j\n",
      "  384.6079+5.8853e-14j  524.8468-5.6523e-14j  524.5176-9.0486e-14j\n",
      "  384.0719-6.0925e-14j  384.6079+1.8811e-14j  265.5464+2.4561e-15j\n",
      "  265.6359-5.0074e-14j  168.8714-9.1821e-15j  168.0086-2.2293e-14j\n",
      "  265.5464+1.5522e-14j  265.6359+5.8174e-15j  168.8714-2.5878e-14j\n",
      "  168.0086-8.6923e-15j   92.671 -4.6094e-14j   93.1292-1.2842e-14j\n",
      "   38.2909-4.0661e-15j   39.6635-1.4006e-14j   -6.1403+1.5693e-14j\n",
      "    7.8542+1.2932e-15j    6.1674+5.0691e-15j   -6.1403-1.2026e-14j\n",
      "    7.8542+6.1556e-15j    6.1674+4.9650e-15j   92.671 -1.6372e-16j\n",
      "   93.1292-1.2402e-14j   38.2909-2.1373e-14j   39.6635+8.6243e-15j]\n"
     ]
    }
   ],
   "source": [
    "print(sq_beigenvals)\n",
    "#we could try to numerically pair them (IF WE have some guarantee that they come in complex conjugate pairs), which I think is only valid\n",
    "# for the TE case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complex eigenvalues of real matrices\n",
    "Statement: eigenvalues of real matrices ALWAYS come in complex conjugate pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
